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Abstract 


The use of the Global Positioning System for position and attitude determination is 
evaluated for an automated rendezvous and docking mission. The typical mission 
scenario involves the chaser docking with the target for resupply or repair purposes, and 
is divided into three sections. During the homing phase, the chaser utilizes Coarse 
Acquisition pseudorange data to approach the target; guidance laws for this stage are 
investigated. In the second phase, differential carrier phase positioning is utilized. The 
chaser must maintain a quasi-constant distance from the target, in order to resolve the 
initial integer ambiguities. Once the ambiguities are determined, the terminal phase is 
entered, and the rendezvous is completed with continuous carrier phase tracking. 

Attitude knowledge is maintained in all phases through the use of the carrier phase 
observable. A Kalman filter is utilized to estimate all states from the noisy measurement 
data. The effects of Selective Availability and cycle slips are also investigated. 
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1 . Introduction 


1.1. The Global Positioning System 

1.1.1. System Overview 

The Global Positioning System (GPS) is a satellite navigation system that was 
first proposed in the early 1970's. Primarily, the system was to provide United States 
military forces with accurate time and position information anywhere on the globe, at any 
time, and in any weather. Since then, the civilian applications of the system have grown 
and will continue to grow beyond anything the original designers intended. As it was 
originally conceived, the system was to consist of a constellation of 24 satellite vehicles 
(SVs), 21 of which were active with 3 on-orbit spares. With the design and testing phases 
complete, the system is currently being deployed by the United States, and at the time of 
this writing a total of 26 operational satellites, of varying ages and capabilities, are in 
orbit. 

Each GPS satellite broadcasts a unique signal from which a user with the proper 
equipment can determine accurate time and position information. At its most basic level, 
these signals consist of the positions of the individual GPS satellites; thus, the user 
position can be determined from a simple triangulation. In most cases, the user desires 
three dimensional position knowledge, so at least three GPS satellite signals must be in 
view to meet this need. However, because the user clock will generally not be 
synchronized with GPS time, this adds an extra unknown so that at least four GPS 
satellites must be in view for position and time knowledge. Of course, more than four 
satellites may be used to provide a more robust solution. 

The requirement that four GPS satellites be in view dictates the number of 
satellites that must be in the constellation. When deployed in six equally spaced 12 hour 
orbits, the minimum number of satellites considered necessary to provide adequate 
coverage is 21. A measure of the coverage of a satellite constellation is its constellation 
value , which represents the fraction of the Earth and time in which at least four satellites 
will be in view; for the GPS 21 satellite primary configuration, this value is not to fall 
below 0.9960. 1 Thus, access to accurate time and position information is available to 
users over the entire globe, virtually 24 hours a day, in any conditions. 

1.1.2. The GPS Signal 

The GPS receiver determines the measured travel time of the signal from the GPS 
satellite through the use of two pseudorandom noise (PRN) codes. The first, the C/A- 
code (Coarse Acquisition-code), is intended for all users; it is also called the Standard 
Positioning Service (SPS). The second, the P-code (Precision-code) is intended only for 
military or other authorized users. The rationale for the two separate codes is security. 
While the effective wavelength of the C/A-code is 300 m, the P-code is 30 m; thus the 
civilian C/A-code performance is significantly degraded from the military P-code. In 
addition, the GPS signal can be intentionally degraded under a policy called Selective 
Availability (SA), at which point only users with authorized codes can eliminate the 
added error. This scheme may be invoked, for example, in time of military action so that 
the enemy may not benefit from the service. 2 
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The GPS signals are broadcast on two carrier frequencies, LI (1575.42 MHz) and 
L2 (1227.60 MHz), using spread spectrum techniques so that they are less susceptible to 
jamming. Both the C/A- and P-codes are modulated on the LI carrier; however, only the 
P-code is on the L2 carrier. Again, the intention is to deliberately degrade the C/A 
performance. Access to two frequencies permits the user to eliminate ionospheric errors, 
so the P-code solution will be more precise. 3 

The information contained in the GPS signal, called the navigation message, 
provides information on the broadcasting satellite, as well as the entire constellation. 

Over half of the navigation message is dedicated to accurately describing the state of the 
satellite that broadcast the signal. This includes information on the health of the satellite, 
the satellite position and velocity, and information to correct for various errors. The 
remainder of the message contains ephemerides for the rest of the GPS satellites (for 
quick visibility checks and signal acquirement), health status of the GPS satellites, as well 
as information reserved for military use. 

The control of the GPS SVs is maintained through a worldwide network of 
ground stations, with the master control station located at the Consolidated Space 
Operations Center (CSOC) at Falcon AFB, Colorado Springs, Colorado. The main 
responsibility of the ground segment is to ensure that the information broadcast by the 
satellites is within operational specifications. Once per day, the monitor stations upload 
ephemerides and clock information, as determined by the master control station. This, 
combined with the fact that the satellite clocks are accurate to a few parts in 10 n over one 
day, ensures that the SVs remain synchronized, and can thus provide accurate positioning 
information. 

1.1.3. The GPS Observables 
1.1. 3.1. Pseudorange 

The pseudorange is a timing measurement related to both the signal propagation 
delay from the satellite to the receiver, and the clock offsets. If the user clock were 
precisely synchronized with GPS time, then the range to the GPS satellite would be 
determined by multiplying the propagation delay time (At) with the speed of light (c) 

R = cAt ( 1 ) 

The delay time is calculated through a code correlation of the GPS signal and the internal 
code generated by the receiver. However, the user and GPS clocks will generally not be 
synchronized, so the correlation time (At) consists of true travel time plus a clock offset 
component 


At = U, rut +K lwk ( 2 ) 

Thus, the measurement consists of the true range plus a component due to the clock 
offsets, hence the name pseudorange 

R = c{At, rm +At rhck ) (3) 
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As mentioned above, there are four quantities of interest, three positions and the clock 
offset. Thus, the four pseudorange equations corresponding to the four GPS satellites 
must be solved simultaneously; this will yield the user position and clock offset. This 
process is explored in more detail in Section 2.6, Position Triangulation. 

1.1. 3. 2. Carrier Phase 

The carrier phase is a second observable available to users of GPS, possessing 
both advantages and disadvantages. On the positive side, because the GPS wavelength is 
approximately 19 cm, this results in a significant increase in accuracy over the use of 
C/A-code or even P-code. However, the implementation of the carrier phase is both more 
complicated and more sensitive. An event called a cycle slip can occur in practice; this 
happens when the receiver miscounts the carrier phase, leading to errors in the position 
knowledge. Still, the general consensus is that the advantages to be gained outweigh the 
added difficulties. 

Consider a GPS carrier phase receiver, immediately after carrier tracking has 
begun. Due to the cyclical nature of the phase, the instantaneous phase measurement will 
be in the range 0 to 271. However, the total phase measurement is given by this value plus 
some integer number of wavelengths representing the remaining distance to the GPS 
spacecraft (N). Thus, the phase measurement is described by 

rue + At clock) + N (4) 

where, as with the pseudorange, the measurement reflects both the true range and a range 
associated with the clock offset. If the integer N were known, then complete position 
knowledge could be obtained through the continuous tracking of the phase. In general, 
GPS receivers track the phase by initializing an internal counter to a preset large negative 
integer. Then, at every positive phase zero crossing, this counter is incremented, yielding 
knowledge of the number of wavelengths passed since the measurement was begun. 
Unfortunately, N, termed the initial integer ambiguity, is generally not known. Thus, 
before the carrier phase information can be of any value, the integer ambiguity must be 
resolved. Several methods are available for performing this task, all of which are 
described in more detail in Section 3.2, Differential GPS and the Integer Ambiguity. 

1.1. 3. 3. Differential GPS 

When dealing with two or more GPS receivers in "close" proximity, it is 
important to realize that many of the errors present in the signal are common to both 
receivers when the same GPS satellites are used for triangulation. The definition of close, 
in this sense, must be evaluated on a case by case basis, dependent on the particular 
requirements in question. One example of a common error is the ionospheric and 
tropospheric delays in the GPS signal. In addition. Selective Availability errors can also 
be considered common to receivers relatively close together. When dealing with such 
errors, higher accuracy results can be obtained by utilizing differential GPS. This mode 
of operation entails subtracting the measurements of one receiver (pseudorange or phase) 
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from another, thereby obtaining a relative measure of position. In this way, common 
errors are eliminated in the differencing process. 

1 .2. Scope of the Investigation 

This study is primarily intended to benefit automated rendezvous missions where 
the chaser is some type of unmanned craft, and the target is cooperative. The term 
cooperative in this sense basically means that the target is broadcasting all necessary GPS 
information to the chaser. Most likely, these types of missions will fall into the resupply 
or service areas. For example, an increasingly possible scenario is the use of a small, 
automated craft for the resupply of the International Space Station. Also, the methods 
explored in research could be applied to the capture and repair of a malfunctioning 
satellite, or the refueling of an operative vehicle. In addition, in the wake of the Galileo 
and Mars Observer failures, it has been proposed to utilize some type of on-orbit 
inspection craft to test such things as antenna deployment; in this way, necessary repairs 
can be made before the craft leaves the Earth. Hence, the subjects explored in this 
investigation are applicable to a wide variety of rendezvous scenarios. 

1.2.1. The Nominal Mission 

For the purposes of this investigation, the typical simulation will begin assuming 
the chaser is on-orbit, trailing the target by 50 km. For clarity, a simplified rendezvous 
sequence is shown in Figure 1.2.1. a below. The first phase of the rendezvous, called the 
homing phase, begins with a maneuver (AV,) that enters the chaser on a trajectory 
towards the target. During this phase, C/A pseudoranges will be used to estimate the 
position of the chaser alone, independent of the target. One or more corrective burns 
(AV 2 ) may be performed to account for trajectory errors. At the end of the first phase, it 
is not intended for the chaser to meet the target; rather, the chaser enters a stationkeeping 
mode a nominal distance behind the target (approximately 2 km). To accomplish this, a 
bum (AV 3 ) is performed to reduce relative motion between chaser and target. 


Target 

Q 


Stationkeeping 

Phase 

/ 

AV*<C>AVa 


AV* 

Terminal 

Phase 



Initial 

Chaser 

Position 

+ 

, AVt 


Homing Phase 


Figure 1.2.1. a: The Nominal Rendezvous Scenario (not to scale) 
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In the second phase of the rendezvous, the stationkeeping phase, differential GPS 
with carrier phase positioning is used to gain enhanced relative position knowledge. As 
mentioned previously, when the carrier phase is used, extra unknowns (the integer 
ambiguities) must be resolved. Thus, this phase is required solely to initialize the carrier 
phase positioning. No maneuvers are needed or wanted; a constant chaser position 
relative to the target is desired as much as possible to hasten the determination of the 
unknowns. 

Once the integer ambiguities are resolved, a maneuver is performed to close the 
remaining 2 km between the spacecraft (AV 4 ). At this point, the chaser is in the terminal 
phase of the rendezvous, and any trajectory errors are eliminated through a final 
corrective burn (AV S ). The continuous tracking of the carrier phase will yield highly 
accurate position and velocity information until the rendezvous is complete (AV 6 ). 

During all phases of the mission, carrier phase data is processed to estimate the 
attitude of the chaser. It is important to note that because of the close proximity of the 
receivers, it is not necessary to resolve any integer ambiguities in the attitude 
determination. No attempt has been made to model an attitude control system. 

1.2.2. Simulation Method 

This simulation utilizes the Kalman filter approach to estimate the state from a 
series of noisy measurements. While the Kalman filter is a linear filter, some processes in 
this study (GPS position triangulation) are nonlinear in nature. In addition, the SA signal 
is non-Gaussian, which complicates the interpretation of the simulation statistics. The 
equations of motion are linearized in all cases for use with this type of filter, however, 
whether this approach is valid is a question to be answered by this study. 

This simulation uses a "truth" and an "estimate" model to investigate the system 
of interest. The truth model is based on the propagation of the state via the linearized 
equations of motion (position and attitude), with state noise added to account for 
mismodeling of the system. The Kalman estimate of the state is derived from the noisy 
measurements (pseudorange and phase). Simulation strategies are discussed in more 
detail in Sections 2.7, The Kalman Filter , 3.3, Kinematic on the Fly GPS and Kalman 
Filtering, 4.3, Kinematic on the Fly GPS and Terminal Filtering , and 5.5, The Attitude 
Kalman Filter. 

One final topic involving the data processing in this investigation should be 
discussed. During the homing phase, the simulation only processes data from the chaser; 
the target position and attitude are assumed known. While this obviously will not be the 
case in a true mission, it is assumed that the target will possess GPS sensors. It then 
follows that the position and attitude knowledge demonstrated on the chaser in this 
simulation will be achievable by the target as well. In this case, an estimate of the 
nominal error (chaser and target) during the homing phase can be obtained by the root- 
sum-square (RSS) of the two components; this will lead to an increase in the error 
demonstrated in this investigation by a factor of -J2 . At worst case, the error would be 
the sum of the two components, leading to a factor of 2 increase in the homing phase 
results contained herein. Once the homing phase is complete, the chaser and target 
switch to differential carrier phase positioning, and pseudorange data is no longer 
processed. In this way, higher accuracy relative state information can be obtained. Thus, 
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because no assumptions are made on the position of the target, the factor of V2 is not 
necessary for the final two phases. 

1.2.3. Investigation Objectives 

The fundamental question this research seeks to answer is whether two 
cooperative spacecraft can rendezvous within a prescribed tolerance using GPS alone, i.e. 
unaided by any other navigation or attitude sensing system. Many specific points to be 
answered fall under this general question, including: 

• Can Kalman filtering yield sufficiently accurate results in all phases of the 

rendezvous (pseudorange and carrier phase), for both position and attitude 
sensing? 

• How will Selective Availability affect the pseudorange results? 

• What type of guidance strategy should be used? 

• Can the integer ambiguities be resolved in the orbital situation? 

• Will carrier phase positioning provide sufficient state knowledge for terminal 
rendezvous? 

• Could a carrier phase cycle slip lead to catastrophic failure? If so, can a cycle 
slip be detected and repaired? 

All of these topics are investigated and discussed in the sections below. 


2. Position Determination from Pseudorange Measurements 

2. 1 . Homing Phase Overview 

This phase of the investigation analyzes the position knowledge of the chaser as it 
homes in on the target, covering the range from AV, to AV 3 in Figure 1 .2. 1 .a. It is 
assumed that the orbital insertion is performed satisfactorily, with the chaser trailing the 
target at the specified distance, typically 50 km. Based on the initial knowledge of the 
chaser's position, a maneuver is performed to initiate the rendezvous. C/A code 
pseudoranges will be used to determine the chaser position until it arrives within about 
two kilometers of the target. At user specified intervals, the chaser receives measurement 
updates perturbed by SA. A Kalman filter is utilized to obtain the optimal estimate of the 
chaser's position. These estimates are used to calculate the necessary corrective burns to 
rendezvous with the target. The system models and results are described below. 

2.2. Orbital Mechanics 
2.2.1. Target and Chaser 

This section highlights the mechanics between the target and the chaser. For 
simplicity, the computations are performed in the target orbital frame. As seen below, 
this frame is formed with the X-axis opposite the target velocity vector, the Y-axis in the 
radial direction, and the Z-axis normal to the orbital plane. 
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Earth 


Figure 2.2.1. a: The Target Orbital Frame 

For the target orbital frame, the linearized equations of relative motion of the 
chaser assuming the target is in a circular orbit are 4 
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y + 2n x — 3n 2 y = 0 (5) 
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Thus, given the position and velocity of the chaser relative to the target, the chaser's 
position at a later time t is given to first order by this equation. For this simulation, the 
time period of interest is in between measurements, or t = At . It should be noted that the 
matrix in the above equation is the state transition matrix used in the Kalman filter, this 
will be explored in more detail in later sections. 

2.2.2. Chaser Maneuvers 

Given an initial position jc 0 , y 0 ,z 0 , it is possible to determine the velocities 
necessary to arrive at the origin in a specified time x. The solution to the linearized 
equations can be inverted to yield the following conditions 6 


*0 sinwx +_y 0 [6A7T sinnx -14(l-cos/ix )] 

3x sin/ix -8(l-cos«x )/ n 
2x 0 (\ -cos nx ) + y 0 (4sin«x -3«x cos«x ) 

3x sinrzx -8(l -cosox )/n K 

Z£ o n 
tan/ix 

When there is a burn (such as at the start of the simulation), the true and estimate velocity 
increments must be determined. Based on the above equations, the velocities to arrive at 
the origin are calculated based on the current estimate of position, because the true 
position is not known. The velocity estimate after the maneuver becomes the nominal 
value calculated from the equations above. The true velocity is set to this nominal 
velocity plus terms accounting for errors in the burn. 

To account for maneuver execution errors, two parameters are used which * 
describe duration and direction errors in the burn. Duration, or magnitude, errors can be 
attributed to burn time errors, incorrect specific impulses, and the like. These errors are 
modeled through the use of the standard deviation of the magnitude error, an input 
parameter assumed to be known at the start of the simulation. Direction errors could be 
due to thruster misalignment, spacecraft attitude errors, and so forth. These effects are 
modeled through the standard deviation of the error in the two directions perpendicular to 
the burn axis. A numerical value is obtained from the root-sum-square (RSS) of two 
components. The first is a floor AV direction error, representing thruster misalignment. 
This is an input parameter, and can be changed to investigate a variety of thruster 
accuracies. For this study, this value is nominally set at 0.01°. The second component is 
obtained from the attitude knowledge at the time of the maneuver. The rationale for this 
is that if the attitude is known only to a certain degree, then the accuracy of the pointing 
of the AV cannot be any better than that. The actual procedure to determine the true 
velocity increments is described in Section 2.8, Maneuvering Errors. Through the use of 
these methods, effects of errors in thrusting can be investigated. 

2.2.3. GPS Satellites 

Whereas the calculations for the target and chaser were performed in the orbital 
frame, the GPS processing was carried out in an inertral frame fixed at the Earth's center. 
The GPS satellite orbits were propagated using universal variables 7 . In order to avoid 
having to recalculate the GPS states every simulation, a pre-processing program was 
written that accepts the orbital elements for the constellation and constructs a table for 
subsequent utilization. In this way, the GPS orbit prediction need only be performed once 
for each constellation, thus building a library of GPS state tables. For example, one data 
file may have a 16 satellite constellation, and another may have 21 satellites. The 
constellation used to obtain the results in this investigation consisted of an optimal 21 
satellite configuration 8 . 

In addition, given the orbital elements of the target, its ephemerides are 
propagated and stored via the pre-processing. Any conversion between inertial and local 
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frames is performed through the use of the direction cosine matrix (DCM). For example, 
the inertial position of the chaser in terms of the target orbital frame is 


chaser _ origin 




( 8 ) 


This equation states that the inertial position of the chaser equals the inertial position of 
the origin of the orbital frame plus the DCM from the orbital to inertial coordinates times 
the chaser position in orbital coordinates. For a circular orbit, the DCM is found from the 
relation 



where r is the position of the body in inertial coordinates. Likewise, when converting 
from inertial to orbital coordinates, the orbital position is found from 

(10) 

where the direction cosine matrices are related by 

CD 

Thus, given any inertial position, these relations permit the conversion into orbital 
coordinates, and vice versa. 

2.3. Line of Sight Determination 

GPS receivers have software that selects the optimal set of GPS satellites with 
which to determine the current position. Before this can be carried out, the line of sight 
to all satellites must be determined for the time periods of interest. In addition, it was 
desired to include the capability to set a mask angle from the surface of the Earth 
(nominally 5°), in order to model the reduced visibility through the atmosphere near the 
horizon. This section highlights the method used to determine receiver-satellite line of 
sight. Once the chaser orbital coordinates are transformed into inertial coordinates, the 
following process is carried out for every SV. 

a) If the component of the SV position in the direction of the chaser is greater 
than the radius of the chaser, then the SV must be in sight, and the line of sight 
process is complete; otherwise, continue to step b. The SV position in the 
direction of the chaser is found from 
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(12) 
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b) Determine whether a line connecting the chaser and the SV positions crosses 
the Earth's surface (assumed to be a sphere). If this is true, the SV is not in 
sight, and the line of sight determination process is ended; otherwise continue to 
step c. The line of sight crosses the Earth's surface if the shortest distance 
between the line and the origin (the center of the earth) is less than the radius of 
the Earth. The distance from the origin to the line 
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where a is the direction vector 

a=[.r 2 -X| z 2 -z,] T 


is given by 9 


d = 
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(13) 


(14) 


(15) 


where u = [-x, -y, -z,] T . 

c) From the known receiver and SV positions, construct the following figure: 


Chaser 



Figure 2.3. 1 .a: Line of Sight with Mask Angle Construction 

where a is the mask angle. Three unit vectors are calculated: first, from the 
chaser to the point of tangency on the Earth (a), second, from the chaser along 
the line formed when the Earth tangent is rotated to the desired mask angle (b), 
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and third, from the chaser to the SV (SV, or SV 2 ). The two cross products 
a xSV and SVxb are evaluated; if they are both the same sign, then the SV is 
not in view due to the mask angle. For example, these cross products are the 
same sign for SV, in Figure 2.3. 1 .a, so it is not is sight due to the mask angle. 
On the other hand, SV 2 will yield cross products of opposite signs, so it is in 
sight. This process must be performed for every SV at every time step to yield 
the required line of sight information. 


2.4. SV Selection 

The minimum number of SVs required to determine the user's coordinates via 
pseudorange data without any prior information is four, to account for the three position 
coordinates and the receiver clock offset from GPS time (see Section 2.6, Position 
Triangulation). The choice of the best four satellites out of the n SVs in view is very 
critical, because it directly relates to the accuracy of the solution. The optimal choice is a 
function of several variables, including geometry, signal strength, health of the SV, and 
so forth. Currently, the simulation only takes into account the geometry of any given 
triangulation; however, it is possible to enhance the programming and include these other 
effects as well. 

This simulation utilizes the Geometric Dilution of Precision (GDOP) as the 
standard of selection of any four satellites. The GDOP is basically a "measure of 
goodness" of the geometry of any particular GPS solution; it can be thought of as the 
inverse of the volume formed by the four SVs and the receiver. Thus, if the SVs are all 
almost directly overhead of the receiver, then the volume would be small, and the GDOP 
high, implying a poor geometry. However, if one of the SVs was overhead, and the other 
three were all at a low inclination to the chaser and equally spaced at 120° apart, the 
volume would be close to the maximum, implying a very low, and very desirable, GDOP. 

The formula for calculation of GDOP is 


GDOP = 


1 

3V 


(16) 


where V is the volume of the tetrahedron formed after connecting the endpoints of four 
unit vectors pointing from the receiver to their respective i ,h satellite, b, is the area of the 
i* lateral face of the endpoint tetrahedron, and v, is the volume of the tetrahedron formed 
by the i* lateral face of the endpoint tetrahedron and the receiver location 10 . Note that the 
area of a triangle with the vertices P,, P 2 , P 3 is given by 1 1 

Yi 

( 17 ) 

and the volume of the tetrahedron with the vertices P,, P 2 , P 3 , P 4 is 12 
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It is possible to evaluate every combination of four of the n SVs in view, and then 
choose the combination that yields the minimum GDOP. The number of evaluations of n 
SVs taken 4 at a time is given by the binomial coefficient 13 


N = 


' Si ) 

U«- 4 ) !4! J 


(19) 


so that as the number SVs in view increases above four, the necessary number of function 
evaluations increases rapidly. To reduce this computational requirement, a sub-optimal 
GDOP selection criterion 14 has been incorporated. The process involves determining the 
three separate GPS satellites that satisfy the following criteria. The vector from the 
receiver to the first SV has the greatest component along the radial direction of the 
receiver; the vector from the receiver to the second SV has the greatest component 
opposite the receiver's velocity vector, and the vector from the receiver to the third SV 
has the greatest component perpendicular to the above two vectors. Finally, for each of 
the remaining satellites, the GDOP is calculated with the above three SVs, and the 
resulting minimum is selected. Fundamentally, this procedure ensures that three of the 
four SVs will be chosen so as to form a nearly orthogonal set of axes, resulting in a large 
enclosed volume, and hence a small GDOP. In practice, this procedure achieves GDOP 
values on average about 1.5 times the minimum. This is illustrated in Figure 2.4. 1. a, 
which shows the minimum and sub-optimal minimum GDOP values for a typical 
rendezvous scenario over one orbit of the target. 
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Figure 2.4. 1. a: Typical Minimum and Sub-Optimal GDOP Values 


The most important aspect of this method, however, is that it requires only n - 2 GDOP 
evaluations. The computational savings are shown in Table 2.4. 1. a. 
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Table 2.4. 1. a: Comparison of Optimal and Sub-Optimal SV Selection Efficiencies 

Number of SVs Optimal GDOP Sub-Optimal GDOP 
in Sight Evaluations Evaluations 

5 5 2 

6 15 3 

7 35 4 

8 70 5 

9 126 6 

10 210 7 

A point to emphasize about this SV selection method is that from one time step to 
the next (At = 3 min), the chosen four satellites rarely remain the same. This is a positive 
trait in the sense that it virtually ensures independent measurements at every pseudorange 
epoch. On the other hand, this method requires that the pseudorange receiver maintain a 
lock on all satellites in view. For illustrative purposes, a typical plot of SV usage over 
one orbit is shown on Figure 2.4.2.a. A dot means that the satellite is in view, and a circle 
means that the satellite is in use. 
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Figure 2.4. 2. a: Typical SV Usage 


2.5. Selective Availability 
2.5.1. Overview 

It was desired to investigate the effects that Selective Availability may have on 
performing an automated rendezvous mission; thus, an SA model is included. Numerous 
studies have simulated SA statistically; for example, some models assume that the 
pseudorange error due to SA is an exponentially correlated random variable with a time 
constant of 15 minutes 15 . Other models have attempted to account for both "orbit SA" 
(perturbed GPS SV ephemerides reported to the receiver), and "dither SA" (false drift of 
the SV clocks) 16 . In this case, both effects are modeled as second order Markov 
processes, with correlation times of 3 and 60 minutes, respectively. This simulation 
utilizes an SA model based on actual GPS data developed by Braasch 17 . The advantage 


20 


that this model has above the statistical approaches above is that it was derived from 
actual GPS data. 

2.5.2. SA Residual Generation and Incorporation 

The procedure for simulating SA is begun by obtaining an actual sample of 
pseudorange data from one GPS satellite. With only one SV, the issue of the receiver 
clock offset from GPS time poses a problem. If a highly stable crystal oscillator is used, 
the clock offset could be represented by a first order polynomial. Then, the residuals are 
calculated by subtracting the best fitting straight line from the collected data. Braasch 
states that these residuals primarily represent SA because the variance due to SA is two 
orders of magnitude greater than all other error sources. It is important to realize that this 
method will not identify a non zero mean component of SA. 

Once the SA residual time series is collected, an optimum filter G(oo,t) can be 
constructed such that when the data is passed through it, the output is white noise with 
variance a 2 . Next, the inverse of G(co,t), H(a),t), is determined. Finally, statistically 
equivalent SA residual data is acquired by passing white noise with variance a 2 through 
the inverse filter. Fundamentally, this procedure reduces to producing a time series 
whose power spectrum is the same as the original signal. In this way, an arbitrary amount 
of SA data can be generated for the purpose of this simulation. In addition, it is very easy 
to modify the filter parameters should new pseudorange data indicate a change in the 
nature of the SA signal. 

The incorporation of this model into the simulation is straightforward. The 
necessary length of SA data for each GPS satellite is calculated according to the 
procedure described above. This simulation investigates two levels of SA measurement 
noise, 27 m and 39 m, as described by Braasch. Then, when the GPS pseudoranges are 
calculated, the SA error corresponding to the correct SV and time is simply added to that 
pseudorange. These perturbed pseudoranges are then entered into the triangulation 
routine to determine the position estimate, which is biased by SA. These are the 
measurements used by the Kalman filter to optimally estimate the state. 

2.5.3. The Non-Stationarity of Selective Availability 

A process is said to be stationary if its statistical properties do not vary with the 
passing of time. The model used in this investigation assumes that the SA residual data is 
stationary. However, in a more recent paper 18 additional SA data samples were analyzed, 
and it was determined that several of the records did exhibit non-station ary 
characteristics. Unfortunately, this complicates the modeling of SA; instead of simply 
modeling a system with constant coefficients, the coefficients now change with time. On 
a positive note, the recently collected data did exhibit stationarity for periods of up to one 
and a half hours. As a result, it will be assumed that the current SA model will be valid 
for the time periods of interest. 

2.6. Position Triangulation 

The final operation involves the position triangulation from GPS pseudoranges; 
this basically follows the standard procedure 19 . With the four GPS satellite positions (X j? 
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Y,, Z,) and the four simultaneous pseudorange measurements (R^, the user position (x, y, 
z) and clock bias (At) can be determined from the simultaneous solution of 


ta-«4/) J =(*-X,)’+(y-i;)‘+U-Z,) ! i— 1...4 (20) 


The solution involves linearizing these equations, and applying Newton's method. Given 
an a priori estimate of the chaser's state ( x,y,z,At ) the corrections to this a priori 
estimate (5x, 8y, 5z, 5At) can be computed by first expanding the measurement equations 
and keeping only first order terms, so that 


_ a 3/? 5 dR ~ 3 R t ~ dR , A 

R ‘ =R ' + li &x+ t y 1?* *£ 


( 21 ) 


The partial derivatives are given by 
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Thus, the problem can be formulated as the matrix equation 
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so that 
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Therefore, the corrections to the a priori estimate 8u can be solved for iteratively, until 
the desired degree of accuracy is achieved. Incidentally, the A-' matrix is a direct 
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indicator of the instantaneous geometry of the system, and will be used in the Kalman 
filter to properly transform the covariance; this is discussed in more detail in the next 
section. 

2.7. The Kalman Filter 
2.7.1. Filtering Procedure 

Knowing that the GPS pseudoranges have been perturbed by SA and other errors, 
it is desired to filter these noisy measurements to obtain an optimal estimate of the 
spacecraft position. A standard Kalman filter is utilized to acquire this estimate, based on 
the development in Gelb 20 . As outlined in the introduction, the Kalman filter utilizes a 
"truth" model and a "knowledge" model. Both models are propagated through the use of 
the state transition matrix, <I>; however, the truth model accounts for errors in system 
modeling by adding a random vector w, with zero mean and covariance Q k (the state 
noise). This covariance matrix is an input parameter, so that a variety of system errors 
can be investigated. The truth model is of the form 

x* = + w*., w, * W(0,QJ (27) 

Note that the state transition matrix, is the same as described in Section 2.2.1, Target 
and Chaser, except that the state is augmented to account for the receiver clock bias. 

This investigation assumes that over the time period of interest, there is no drift of the 
user clocks. Thus, the clock bias term in the state transition matrix is unity. 

The measurement for this model is the receiver position calculated from the GPS 
signals. The measurement equation is of the form 

z*=r + v* y k «Af(0,Rj (28) 

It should be made clear that r in the above equation is the position estimate perturbed by 
SA only. This is calculated by taking the true pseudoranges from the receiver to the GPS 
satellites and adding the Selective Availability error. Then, these pseudoranges are used 
to determine the perturbed position measurement, based on the procedure outlined in 
Section 2.6, Position Triangulation\ this procedure yields the SA perturbed position 
vector r . This model also includes a noise term v with mean zero and covariance R k (the 
measurement noise), yielding the noisy measurement vector z. As with the state noise 
matrix (Q k ), R k is an input parameter to allow for flexibility in modeling the system. 

The measurement noise matrix R k must be expressed in terms of the pseudorange 
errors, because the pseudorange error lies along the line connecting the receiver and the 
GPS satellite. The geometry of each individual measurement is then converted to inertial 
coordinates using the A-' matrix described in the previous section. The proper R k is 
calculated from 


R*= A 'R t A“ ,t = g 2 A~'I A“‘ t = <7 2 A -i A~ |T (29) 
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where R t has been assumed to be diagonal (i.e., spherical in pseudorange), with 

Xpseudomnge 

noise variance levels of a 2 . 

As with the truth model, the estimate model utilizes the state transition matrix to 
propagate the estimate; however, the procedure is slightly different. Starting after a 
measurement update, the state is propagated until the next measurement via the relation 

< 30 > 

where (+) corresponds to the time after a measurement, and (-) is the time immediately 
before the next measurement. The error covariance is projected in between 
measurements by 


< 31 > 

Upon the reception of another measurement, the state and covariance estimates are 
updated to account for this new information. First, the Kalman gain matrix is computed 

K, =P,(-)H;[H,P,(-)Hl +R,] H (32) 

Then, the state and covariance updates are 

i,(+) = i»(-)+K.[ I *- H * i *(-)l (33) 

P 1 (+) = [I-K,H,]P,(-) (34) 

where these new states reflect the new information contained in the latest measurements. 
This entire procedure is repeated until the specified end time. 

2.7.2. Starting Requirements 

Before the simulation process can begin, there are several input parameters which 
must be specified. This section lists these parameters, and also includes some typical 
ranges of values, in order to provide a more clear indication of the nature of the analysis. 
To begin, an initial nominal state (x 0 ) must be specified, as well as an initial state 
covariance matrix (P 0 ). The chaser was chosen to be trailing the target at the start of the 
simulation by nominally 50 km; in the homing phase, the chaser is to proceed on a 
trajectory to 2 km in front or behind the target (an input parameter), at which point the 
differential GPS processing will begin. P 0 is assumed to be diagonal, with position and 
velocity standard deviations of 60 m and 1 m/s, respectively. Assuming that the chaser 
has been tracking GPS since orbit injection, the position knowledge should be around 30 
m RSS; a worst case value of twice this was taken for this simulation. From this nominal 
initial state, values for the true and estimated states are determined by calculating two 
different error vectors consistent with P 0 , and adding them to the nominal state. 
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In order to propagate the true state, the state noise covariance matrix (Q k ) must be 
specified. This matrix is assumed to be diagonal, and to be a fraction of the drag force at 
the rendezvous altitude. The standard deviation of the diagonal terms is calculated by 
first determining the acceleration due to drag at the spacecraft altitude 21 
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D ~ 


P 

2 



(35) 


where p is the density at altitude, and V is the velocity of the spacecraft. The reciprocal 
of the terms in the parenthesis is called the ballistic coefficient, and is an input parameter 
for the simulation. For this investigation a value of 90.9 kg/m 2 (18.6 lb/ft 2 ) is assumed; 
this is calculated using a mass of 200 kg, a cross sectional area of 1 m 2 , and a drag 
coefficient of 2.2. This choice reflects the fact that the nature of the chaser is a small 
resupply or repair spacecraft. The above acceleration is integrated once over the specified 
At to obtain an "uncertainty velocity due to drag", and twice to obtain an "uncertainty 
position due to drag". These values are then used as the standard deviation of the position 
and velocity in the state noise matrix 
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where f r and f v are scale factors, nominally set to 0. 1 to represent typical errors in pC D . 
This procedure yields the state noise matrix for this simulation, and will serve to reflect 
uncertainties in modeling the system. 

The last parameter which must be specified is the measurement noise. As was 
described in the above section, the measurement noise matrix must be expressed in terms 
of the pseudorange measurement errors. This is because these errors lie on the path 
between the receiver and the GPS satellite. The standard deviation of these pseudorange 
errors are assumed to be 27 m and 39 m for this study. As described previously, these 
errors are then transformed through the use of the A-' matrix for each particular 
measurement geometry. 

2.7.3. Correlated Pseudorange Measurement Errors 

An assumption in the Kalman filter is that the measurement errors are 
uncorrelated. This assumption places a restriction on the time permitted between 
pseudorange observations. More clearly, if pseudorange measurements are taken at one 
second intervals during the simulation, it is assumed that these would be correlated. 
Because a fundamental assumption is violated in this case, the filter will not be an 
unbiased, minimum variance, consistent estimator 22 . In order to avoid this problem, the 
time in between measurements is nominally set to three minutes (private communication, 
Kathy Thornton, Jet Propulsion Laboratory, September, 1993). Over this span of time, 
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the assumption is that the observations are not correlated, preserving the optimality of the 
Kalman filter. 

2.7.4. Covariance Consistent, Random Vectors 

One last point must be mentioned about the use of the Kalman filter. The method 
at various points requires the generation of error vectors consistent with a given 
covariance matrix. For example, state and measurement noise vectors must be generated 
such that they are consistent with the matrices Q k and R k , respectively. To generate such 
random vectors, an eigen-analysis is performed on the covariance matrix in question. The 
square roots of the eigenvalues (the standard deviations) are each multiplied by a different 
normally distributed random number. These products are then multiplied by the 
corresponding eigenvectors, and the resulting vectors are added to yield the statistically 
consistent random vector. This is expressed mathematically as 
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where o- are the square roots of the eigenvalues, ^ are normally distributed random 
numbers with mean zero and unity variance, and v jk are the components of the 
eigenvectors. 


2.8. Maneuvering Errors 

In order to investigate possible control laws for the rendezvous maneuver, the 
capability of modeling spacecraft maneuvers must be incorporated in the simulation. 

This was described in Section 2.2.2, Chaser Maneuvers. In addition, the ability to 
simulate errors in these corrective bums has been included as well. The modeling of bum 
errors is performed through the use of two variables that account for magnitude and 
direction errors. As mentioned previously, magnitude errors could be due to burn time 
errors, specific impulse errors, and the like, whereas di rection errors could be attributed to 
spacecraft attitude errors, thruster misalignment, or other effects. Specifically, the 
quantities chosen to model these burn errors are the standard deviations of the magnitude 
and direction of the burn as a percentage of the burn itself. Thus, as one would expect, 
with a larger AV, there will be associated a larger burn error. These variables are input 
parameters, so that a variety of configurations may be investigated. 

For this research it is assumed that all maneuvers are instantaneous. This implies 
that the position and its covariance do not change during the maneuver, whereas the 
velocity and its covariance do change. The actual procedure to determine the new 
velocity values involves the following steps. The desired velocity, V^, is determined via 
the method described in 2.2.2, Chaser Maneuvers , and the nominal AV is calculated 


A V = V - V 
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The velocity estimate after the burn becomes this desired velocity, because lacking any 
other knowledge, the new velocity should on average be the nominal expected velocity 

V = V*, (39) 

Next, the corrupted truth velocity after the maneuver is calculated. As mentioned 
above, there are two parameters describing the AV error. The first, a m , is the standard 
deviation of the magnitude of the corrective burn, given as a percent of nominal standard 
deviation. Thus, the component of the perturbed AV vector along the nominal AV axis is 
given by 


AVHlAVjl + 'ia.) (40) 

where r, is a normally distributed random variable with zero mean and unity variance. In 
this study, a value of 0.5% is utilized for the nominal a m . 

The second parameter, a d , describes the direction deviation of the AV vector from 
the nominal. It represents the standard deviation of the bum error along two orthogonal 
axis, each perpendicular to AV,. As described previously, the numerical value for CT d is 
obtained from the RSS of two components. The first component accounts for thruster 
misalignment, and can be considered as a floor value for the direction error; this is 
nominally set at 0.01°. The second contribution is derived from the knowledge of the 
chaser's attitude (Section 5, GPS Attitude Determination). The square root of the trace of 
the attitude covariance matrix is taken as the angular measurement error. To convert this 
to a percent of nominal AV, use the relation 

OdU* = tan ( -\/Tr ace ( P aUiUuk ) ) • 1 00 (41) 

Then the AV components in the cross directions are determined from 

A'HlAVjw,) 

AVHlAVjr.o,) 


Thus, these three components form a perturbed AV vector in a frame attached to the 
nominal AV vector 


AV av = AV„, 
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T 2®d 


(43) 


This vector represents the true velocity after the maneuver, expressed in a coordinate 
frame attached to the nominal velocity vector. All that remains is to convert this vector to 
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the target orbital frame, utilizing standard transformation procedures. Finally, the true 
velocity is found from the relation 




( 44 ) 


where cp AV ^ (f)/ is the DCM from the "nominal velocity frame" to the target orbital frame. 
To summarize, the new velocity estimate after a corrective maneuver is the desired value 
of the velocity to arrive at the target in the specified time. The new true velocity is this 
desired value plus extra terms accounting for possible execution errors. 

Along with an instantaneous change in the velocity vector, there must also be an 
instantaneous change in the velocity covariance. Knowing that the parameters a m and a d 
describe the burn errors in the nominal AV frame, the covariance matrix in this frame is 
given by 


L AV 


(lAVj|c„) ! 

0 

0 


0 

flAVjK,)’ 


0 


0 

0 
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(45) 


The conversion of this covariance to the target orbital frame is carried out through the use 
of the covariance propagation law 

of ={ Pa V-Mo/ ^*AV AV (46) 

Finally, the velocity covariance after a corrective maneuver is given by adding this new 
contribution to the old value 


- p +p 

r „ld ^ 1 lof 


(47) 


2.9. Pseudorange Results 
2.9. 1 . Simulation Integrity 

Before proceeding to the full scale analysis, it was desired to investigate some of 
the fundamental outputs to validate the simulation. Several types of data were examined, 
all of which were verified through intuitive analyses. This was performed simply to 
achieve order of magnitude agreement. Shown below are the results from two such 
investigations. It is expected that a smaller GPS constellation will yield fewer SVs in 
sight at any given time; this obvious result is verified in Figure 2.9.1. a. Also, the 
imposition of an elevation angle implies that fewer SVs will be in sight, as shown in 
Figure 2.9. l.b. 
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Figure 2.9.1. a: Effect of Constellation Size on SV Observability 



Figure 2.9. l.b: Effect of Elevation Angle on SV Observability 

2.9.2. Position Knowledge and Monte Carlo Verification 

For linear problems, Kalman filters produce exact results that, if programmed 
correctly, do not require verification. However, some procedures in this analysis are 
nonlinear, so there is a question as to the validity of a linear covariance analysis. As a 
result, it was desired to utilize a Monte Carlo analysis to examine the performance of this 
simulation. Table 2.9.2.a shows the results of four Monte Carlo simulations with 500 
runs each, and their corresponding error analyses. The first column is the standard 
deviation of the measurement noise for the simulations. The remaining two columns 
represent the true position errors as determined from the Monte Carlo simulation, and the 
knowledge of the position errors calculated from the error analysis, both at the end of the 
homing phase. 

Table 2.9.2.a: Monte Carlo Verification of Position 

1 ct Noise (m) Monte Carlo (m) Error Analysis (m) 


NoSA 

15.0 

9.7 

10.7 

Gaussian SA 

27.0 

17.1 

18.3 

SA via Braasch Model 

27.0 

21.6 

19.1 

SA via Braasch Model 

38.3 

28.2 

26.4 


29 





The first case, no SA, represents when the intentional degradation of the GPS 
signal is not activated, and typical measurement errors are assumed to be normally 
distributed with a standard deviation of 15 m. The second case, Gaussian SA, is a test 
case in which the SA is active, and is assumed to be normally distributed with a standard 
deviation of 27 m. The remaining two rows assume SA is activated and consistent with 
the Braasch model, with measurement noise levels of 27.0 and 38.3 m. 

The table shows that for the no SA and Gaussian SA cases, the Monte Carlo and 
error analysis statistics are consistent as expected; these two cases serve to verify the use 
and operation of the linear filter. In addition, the cases involving SA residuals derived 
from real data are also statistically consistent, implying that the assumptions in this 
research (particularly the measurement time interval of 3 minutes) are proper. Thus, the 
Kalman filter can be utilized to estimate the chaser state, even in the presence of 
intentional degradation of the GPS signal. 

2.9.3. Homing Phase State Knowledge, SA Active 

Before presenting the results of this analysis, a brief explanation is in order about 
the presentation of the simulation results in this and the remaining sections of the paper. 

In most cases, results are desired comparing the receiver knowledge of the state to the 
actual state errors. As a result, in most figures, two sets of data are plotted. First, the 
standard deviation calculated in the error analysis is plotted using solid lines; this 
represents the state knowledge. Second, the residual, or the true value minus the 
estimate, is plotted using dashed lines. Also, in figures showing the position and 
velocity, the three components of each are all shown. If all the assumptions in the 
simulation are correct, the error analysis results represent the average of an infinite 
number of simulation residuals. Thus, it can be expected that the error analysis and each 
particular instance of simulation residuals will be consistent. For example, consider the 
rendezvous case with corrective burns at 15 and 84 minutes, and assume SA noise levels 
of 27 m; the time histories of position and velocity knowledge are shown in the figures 
below. In Figure 2.9.3.a, the position knowledge (three solid lines) and the position 
residual (three dashed lines) remain consistent, displaying the same trend. This is the 
case in the velocity knowledge as well, shown in Figure 2.9.3.b. Hence, the simulation is 
statistically consistent, and the Kalman filter can be used to estimate the state in the 
homing phase of an automated rendezvous. For SA noise levels of 27 m, position 
knowledge is maintained on the order of 20 m, and velocity on the order of 2T0 2 m/s 
(RSS), and for SA noise of 39 m, 27 m position and 4- 10 2 m/s velocity RSS knowledge is 
achievable. 
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Figure 2. 9. 3. a: Homing Phase Position Knowledge and Residual 



Simulation Time (min) 

Figure 2.9.3.b: Homing Phase Velocity Knowledge and Residual 
2.9.4. Guidance Strategies 

When evaluating guidance strategies for the homing phase, the final solution will 
be a compromise between two conflicting requirements. First, as little fuel as possible 
(within requirements) should be utilized, and second, the deviation between nominal and 
true final positions should meet the prescribed requirements. It might be argued that the 
second requirement is not critical, because no matter how poor the guidance is during the 
homing phase, good control during the stationkeeping and terminal phases will result in a 
satisfactory rendezvous. However, the flaw in this argument involves the possibility of 
collision. Because the two spacecraft are in such close proximity at the end of the 
homing phase (2 km), significant maneuvering error could lead to mission failure. As a 
result, maneuvering errors must be sufficiently small that the probability of collision is 
within prescribed limits. 

Several simulations with SA noise levels of 27 m were performed to generate an 
example of these concepts. The results are shown in Figure 2.9.4. a, where exaggerated 
system errors were used to illustrate the effects more clearly. It is a plot of the AV 
requirement and the corresponding control errors for varying maneuver times. The 
propellant requirement is expressed as the "statistical AV", calculated from the relation 
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AV„„,,„, = AV, + AV )t , - AV„ — (48) 

where AV^ + AV 3o is the mean plus 3a value of the AV, calculated from Monte Carlo 

simulations, and AV^^, is the AV requirement along the nominal trajectory, in the 
absence of any errors. For the case of the chaser trailing the target by 50 km, the nominal 
AV requirement is 6.03 m/s. The statistical AV represents the excess fuel required above 
the nominal to account for 3a errors. The control errors are deviations in the final 
position of the spacecraft due to errors in maneuvering. 



Figure 2.9.4.a: Exaggerated One Burn AV Requirements and Control Errors 

Several facts can be gathered from this figure. Early maneuvers require smaller 
AV's; however, poorer state knowledge at this earlier time (due to fewer measurements) 
leads to execution errors. In addition, errors in an early maneuver map into larger errors. 
On the other hand, if the maneuver is performed later, the AV requirement increases as 
the chaser approaches the target. There are advantages to a late maneuver, however, one 
being that because the state knowledge has improved (due to more measurements), the 
burn is less likely to be incorrect. Also, a late maneuver implies that any AV error will 
not have the time to propagate into significant control errors. 

It is evident from this analysis that one maneuver is not sufficient to achieve both 
minimum fuel and minimum control error. As a result, several combinations of two 
corrective maneuvers have been investigated. A maximum of two maneuvers was 
utilized in order to maintain a relatively low mission complexity. A good compromise 
seems to be achieved when there is a burn shortly after the homing phase initiation, and a 
bum shortly before the end of the homing phase. This makes intuitive sense based on the 
following argument. Because the state knowledge at the start of the simulation is poor, 
the initial bum, AV,, is in error. After only a few minutes of tracking, the state 
knowledge is improved, so the majority of the control errors of AV, can be eliminated by 
a first corrective burn, AV 2a . An added benefit is that this early burn has a lower 
propellant requirement. Then, with the chaser much closer to the nominal trajectory, a 
maneuver near the final time (AV, b ) is used to ensure that the control errors are within the 
desired limits. 
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The tables below show the results of the investigation into the use of several 
possible two maneuver combinations, assuming SA noise levels of 27 m. Table 2.9.4.a 
investigates varying first corrective maneuver times (AV 2a ), and Table 2.9.4.b varying 
second corrective maneuver times (AV 2b ). Both tables include a AV analysis and a final 
position analysis. As described above, the "statistical AV" is used as the performance 
index for the AV analysis. For the control analysis, the concept of the b-plane is used. 
The b-plane is an imaginary plane perpendicular to the nominal trajectory at the desired 
end time (the end of the homing phase). Thus, any chaser position errors perpendicular 
to the b-plane are time of arrival errors, and do not contribute to a miss of the desired 
terminal point. However, errors in the b-plane do represent rendezvous trajectory errors. 
Therefore, the parameter used for the control analysis is the RSS of the position errors in 
the b-plane at the rendezvous time; this measure describes how far the chaser is likely to 
miss the desired end point. 

Table 2. 9. 4. a shows the effects of varying the first corrective maneuver time, 
while keeping the final burn time fixed at 81 minutes. For this investigation, the total 
homing phase time is 90 minutes. There is a minimum AV requirement at a first burn 
time around 15-21 minutes. Note that, because the final burn time is fixed, this 
establishes the control errors at the homing phase end to around 30 m for all cases, so this 
is not a factor in these results. 

Table 2. 9. 4. a: Homing Phase Two Burn Combinations: First Burn Times 


AV Analysis (m/s) B-Plane Control (m) 


Burn Times 

Mean 

la 

AV 0 + AV, - AV 

Mean 

la 

12, 81 min 

13.0 

4.35 

20.0 

3.16 

29.1 

15, 81 min 

11.2 

3.03 

14.3 

1.93 

30.0 

21,81 min 

11.1 

3.08 

14.3 

1.58 

29.7 

30, 81 min 

11.1 

3.58 

15.8 

0.847 

34.7 

39, 81 min 

13.1 

4.89 

21.7 

2.79 

29.8 

45, 81 min 

15.0 

6.58 

28.3 

2.74 

29.6 

54, 81 min 

18.2 

9.72 

41.3 

4.25 

29.2 


It is also desired to investigate varying the final burn time while keeping the initial 
maneuver time fixed. The initial burn time was chosen to be 15 minutes based on the 
results of Table 2. 9. 4. a. The results, shown in Table 2.9.4.b, indicate that there is a 
minimum AV requirement at a second burn time of 72 minutes. However, unlike above, 
the control plays an important role in this analysis, because the time for minimum 
maneuver requirement does not correspond to the time for minimum final control. A 
typical tradeoff is evident; the final decision depends on whether minimum propellant 
mass or minimum control error is more critical. A reasonable compromise is reached if 
the first corrective maneuver is performed in the neighborhood of 15-20 minutes, and the 
second around 75-85 minutes. Given a more well defined system, including a cost 
function involving velocity increments (AV) and control errors (a), optimization methods 
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could be utilized to determine the minimum propellant and minimum control solutions 
subject to the system constraints. 


Table 2.9.4.b: Homing Phase Two Bum Combinations: Second Bum Times 


AY Analysis (m/s) B-Plane Control (m) 


Bum Times 

Mean 

la 

AV 0 + AV 3o - AV^ . 

Mean 

la 

15, 87 min 

16.1 

7.14 

31.48 

3.64 

23.0 

15, 84 min 

12.9 

4.46 

20.3 

5.68 

26.1 

15, 78 min 

10.2 

2.35 

11.2 

2.11 

40.9 

15, 72 min 

9.37 

2.06 

9.52 

1.31 

49.3 

15, 63 min 

9.41 

2.25 

10.13 

3.89 

76.7 


2.9.5. Effects of Modeling SA Measurement Noise Incorrectly 

For all of the above analyses, the measurement noise statistics were assumed to be 
consistent with the noise levels due to SA. In other words, if the actual standard 
deviation of the SA residual data was 27 m, the measurement noise was also set to this 
level. However, it may not always be possible to estimate the SA noise accurately; as a 
result, it was desired to investigate errors in modeling this parameter. To accomplish this, 
the true SA noise level was fixed at 27 m, and the filter measurement noise was varied 
from 7 m to 47 m. The results are shown in Table 2.9.5.a, which shows both the Monte 
Carlo final position standard deviation and the error analysis position knowledge as a 
function of the assumed measurement noise levels. As discussed in Section 2.9.2, 
Position Knowledge and Monte Ccirlo Verification, when the true and filter statistics are 
the same (27 m), the Monte Carlo and error analysis results are consistent. In addition, 
because the true noise levels on all the simulations are the same (27 m), the Monte Carlo 
results are all on the same level (~ 20 m). On the other hand, as the assumed 
measurement noise is varied, the error analysis values change directly; if the SA 
measurement noise is too low, then the state errors are underestimated, and if the SA 
noise is too high, then the errors are overestimated. For example, if the SA measurement 
noise is assumed to be 7 m, then the error analysis will imply that the position is known 
to within 5.3 m, when in actuality, true errors are on the order of 21 m. This is 
undesirable because unknown errors lead directly to an increase in mission risk. On the 
other hand, if the SA noise is to be 47 m, then the error analysis will yield position 
knowledge of 29.8 m, when the true errors are actually much lower. As a result, the 
correct SA noise level must be utilized if the state knowledge is to be consistent with true 
errors; unfortunately, the nature of the SA noise level is never known ahead of time. One 
possible solution may be to deliberately overestimate typical SA noise levels; by doing 
so, the chances of underestimating the position errors is reduced. 


34 


Table 2. 9. 5. a: SA Measurement Noise Analysis 


Assumed 
Measurement 
Noise (m) 

7 

17 

27 

37 

47 


Monte Carlo 

(m) 

21.4 

19.7 

21.6 

19.8 
19.3 


Error Analysis 

(m) 

5.3 

11.0 

19.1 

23.7 

29.8 


2.10. Homing Phase Conclusions 

This portion of the investigation has brought to light several facets of utilizing 
GPS for the homing phase of an automated rendezvous mission. First, the use and 
operation of the Kalman filter was verified. Test cases assuming Gaussian SA noise were 
shown to be statistically consistent, as were cases utilizing modeled SA noise. Second, in 
the presence of 27 m SA noise, the Kalman filter was shown to yield position knowledge 
in the homing phase on the order of 20 m RSS, and velocity knowledge on the order of 
21 0- 2 m/s RSS. When 39 m SA noise is assumed, the position and velocity knowledge 
approaches 27 m and 4 TO- 2 m/s RSS, respectively. In later sections, it will be shown that 
the stationkeeping mode requires only about 10 minutes, so during this time period, the 
residual velocity maps into a displacement of only about 15-30 m. Third, guidance 
strategies involving one and two corrective maneuvers were investigated. It was found 
that a corrective maneuver performed shortly after the homing phase initiation (-15 
minutes) eliminates much of the error in the initial bum, and requires lower levels of 
propellant than a later maneuver. A second maneuver shortly before the end of the 
homing phase (-10 minutes) ensures that the final control errors will be on the order of 
30 m. Finally, the mismodeling of SA error statistics was shown to lead to 
inconsistencies in the statistical output of the simulation; overestimation of this parameter 
may be advisable to reduce mission risk. 


3. Stationkeeping via Differential GPS 
3.1. Overview 

Section 2, Position Determination from Pseudorange Measurements, described 
the homing phase of the automated rendezvous. During this section of the mission, it was 
only desired to maintain a coarse estimate of the chaser's state, and the analysis showed 
that position and velocity knowledge of around 20 m and 2-10- 2 RSS, respectively, could 
be achieved. Once the chaser reaches about 2 km from the target, a maneuver is 
performed to enter a stationkeeping mode. As with all of the AV maneuvers, errors in the 
burn were modeled as described in Section 2.8, Maneuvering Errors. As a result, the 
chaser does not make a perfect burn to begin stationkeeping, and there remains some 
relative motion between the two spacecraft during this phase. 
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At the outset of the stationkeeping phase, the two spacecraft switch to differential 
GPS using the carrier phase observable to obtain higher accuracy relative positioning; 
pseudorange data is no longer processed. As described previously, when using carrier 
phase positioning, the initial integer ambiguity must be resolved before the phase 
information is of any value. This explains the need for the chaser to establish a quasi- 
constant distance from the target: to resolve the integer ambiguities. If, rather then 
entering a stationkeeping mode, the chaser attempts to resolve the integer ambiguities 
while on a trajectory to rendezvous with the chaser, there is a significant chance that the 
ambiguities will not be determined before the rendezvous time. Without the accurate 
phase data, the chaser will have to rely on the coarse pseudorange data, whereupon the 
risk of mission failure is greatly increased. Thus, the chaser must maintain a quasi- 
constant distance while the GPS carrier is tracked. Filtering the phase measurements, as 
described below, should resolve the unknowns, so that continuous carrier phase tracking 
can be used in the terminal phase of the rendezvous. 

3.2. Differential GPS and the Integer Ambiguity 

This entire development of differential GPS and resolving the integer ambiguities 
on the fly is based on the discussion in Hwang 23 . When utilizing GPS carrier phase 
measurements to determine relative positions, the periodic nature of the signal introduces 
an ambiguity which must be resolved before the entire solution can be determined. 
Consider Figure 3.2. 1 .a, which is an example of a single one dimensional static phase 
measurement situation. 



Figure 3.2.1. a. One Dimensional Static Phase Measurement 24 

The phase difference between points A and B is utilized for the relative positioning. The 
total phase difference is calculated from 


A(j) = (j)^ - <j) B = cos0 • Ax 


(49) 
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This measurement is composed of a part due to the number of whole phase cycles, AiJ)^, 
and a part due to the measured phase difference, A<)) 0 , so that 

A<J) = A<J) 0 + Aij)^ = cos0 • Ax (50) 

The whole cycle count, A$ N , is unknown at the start of the tracking and is called the 
initial integer ambiguity. Thus, the measured phase difference is 

A<t> 0 = cos0 • Ajc- A<|> w (51) 

Only after this ambiguity has been determined does the continuous tracking of the carrier 
signal yield complete knowledge of the relative position of the two receivers. 

It is important to mention that while the above discussion was based on a simple 
single phase difference measurement, in actuality there is a problem associated with this 
type of measurement in the form of the clock error. The true phase difference is given by 

A<J> 0 = cos0 • Ax- Aty N + <oAt (52) 

where the clock error At is due to the time difference in the receiver clocks at A and B. 

To eliminate the clock offset, the phase difference from two different satellites is 
subtracted, forming the double difference measurement. The critical assumption is that 
the measurements are taken simultaneously. Given the following two single difference 
measurements, where the parenthetical superscript denotes the SV, 

A<b 0 (l> =cos0 (1) ■ Ax- A6., (l) +toAf 

(53) 

A<J) 0 <2) =cos0 (2) • Ajc- Ac)),/ 2 ’ +coAf 

the following double difference is formed 

A<|) 0 (,) -A<|) 0 (2) =(cos0 (,) -cos0 (2) )-Ajc + (A<|) w ( 2) -A(J) jv (I) ) (54) 

where the clock offset is seen to cancel. Thus, the measurements used in the simulation 
are the double differences, and the unknown ambiguities are actually the differences of 
two large integers. 

Determining the initial integer ambiguities (actually integer difference 
ambiguities) can be accomplished on the ground typically through three methods. First, 
the ambiguities could be initialized by calibrating the receivers with a known baseline; 
this method provides an instantaneous solution. Second, the antenna swap method could 
be used, in which the receivers exchange positions to resolve the unknowns 25 . Third, a 
static survey could be performed, which would provide a known baseline for future 
surveys. Unfortunately, none of these methods are applicable to the orbital rendezvous 
problem. 

One way to determine the integer ambiguities on-orbit is through the use of some 
type of filtering procedure. There are two main families of filters that could be used, 
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either batch or sequential (the Kalman filter, for this investigation). Noting that the 
homing, terminal, and attitude sections of this study all utilize the Kalman filter (because 
the batch filter is not appropriate for those cases), the same filter was chosen for use in 
the stationkeeping phase, for the ease of implementation. The main disadvantage of this 
approach is that values that are known to be integers are estimated as continuous random 
variables. At the time of this writing, the author was aware of only one investigation into 
the use of a sequential filter that takes advantage of the knowledge that these ambiguities 
are integers 26 . It is possible to utilize a parallel bank of Kalman filters, with each 
corresponding to an integer ambiguity assumption. The filter that corresponds to the 
correct solution converges to 1, while the others converge to 0. A drawback to this 
method is its large computational requirement. For example, real life position uncertainty 
values may be on the order of a few meters, which corresponds to an integer ambiguity 
uncertainty of 100 units in each of three dimensions. The utilization of this method for 
this case would require 1 00-' or a million parallel filters. At this time, this type ot 
computational requirement precludes the possibility of using this approach in practice. 
Thus, a standard Kalman filtering procedure is utilized, as described in the section below. 

It should be noted that the batch filter does possess a particular advantage only for 
the area of resolving integer ambiguities in the stationkeeping mode. If a batch filter were 
used, the chaser would take measurements for some specified time, and then would 
process the entire bank of data, perhaps using a least squares method. The advantage of 
this type of method is that it is easier to utilize the fact that the ambiguities are integers. 
For example, if the least squares analyses yielded non integers for the ambiguities, then 
all possible combinations of integers bounding these values could be investigated. The 
one that provided the minimum residual would be chosen as the solution. The 
implementation of a batch filter for resolving the integer ambiguities in the stationkeeping 
mode requires further study. 

3.3. Kinematic on the Fly GPS and Kalman Filtering 

The idea behind the use of Kinematic on the Fly GPS is basically to utilize a 
Kalman filter to estimate not only the receiver position, but also the initial integer 
ambiguities from the phase double differences 27 . There are two issues that must be 
discussed when utilizing this method. First, Kinematic on the Fly.GPS should not be 
utilized to resolve the integer ambiguities while on a trajectory to rendezvous with the 
chaser (i.e., in the terminal phase); to do so may result in catastrophic failure. If the 
unknowns cannot be determined, then accurate phase position knowledge is not available, 
and only C/A pseudorange data can be used. Large position uncertainty (on the order of 
C/A errors) near rendezvous could result in target-chaser collision. As a result, the 
concept of the stationkeeping phase is utilized; about 2 km from rendezvous, a maneuver 
places the chaser in a quasi-stationary position relative to the target. Phase measurements 
are processed until the ambiguities are resolved and the relative state is fully known. At 
this point, high accuracy differential phase positioning can be used for the terminal phase 
of the rendezvous. 

Second, the observability of the systems in question must be investigated. 
Consider the case where there is no knowledge of the trajectory of the receivers, as in 
Hwang 28 . In this instance, there must eventually be enough measurements to solve for all 


38 


of the unknowns over the time period of interest. For example, consider the case of 
tracking seven satellites (from which six independent double differences can be formed at 
a time), and consider an interest only in the position and the integer ambiguities, and not 
the velocity. For one epoch, or time of measurement, there are three positions and six 
integer difference ambiguities, for a total of nine unknowns. Six measurements cannot 
fully describe these unknowns. However, for two epochs (t 0 and t 0 +At), there are twelve 
measurements, and there are also twelve unknowns: three original positions at t 0 , three 
positions at t 0 +At, and the same six integer ambiguities. Hence, the seven satellite case is 
a well posed problem for position and integer ambiguity knowledge over two epochs. 
Similarly, there are six and five satellite models for these conditions that are both well 
posed; however, the number of epochs required increases above the two needed for the 
seven satellite case. This translates directly into slower filter convergence. On the other 
hand, utilizing larger numbers of SVs becomes impractical because the probability of 
maintaining more than seven in view is reduced. The above discussion can be visualized 
more clearly by inspecting Table 3. 3.1. a. 

Table 3.3. 1. a: Satellite Models and Their Observability, Estimating 
Position and Integer Ambiguity Only 


of SVs 

Epoch # 

# of Msmts. 

# of Unknowns 

5 

1 

4 

9 


2 

8 

12 


3 

12 

15 


4 

16 

18 


5 

20 

21 


6 

24 

24 

6 

1 

5 

9 


2 

10 

12 


3 

15 

15 

7 

1 

6 

9 


2 

12 

12 


For the case of the orbital rendezvous, it is not sufficient to estimate the position 
and integer ambiguities; the velocity must be known as well. This adds three extra 
unknowns at every epoch, so that seven satellites will never be able to sufficiently 
determine the unknowns, no matter how many epochs are used. However, it is still 
possible to utilize the seven satellite model; the key lies in exploiting the information 
inherent in the system dynamics. The assumption of not knowing the trajectory of the 
receiver is valid in most cases; for example, if a GPS receiver were installed on a car or 
an aircraft, the dynamics would be controlled solely by the pilot, and hence could not be 
modeled for use in a Kalman filter. However, in the case of satellites in orbit, the 
trajectories are well known and can be quantified in the state transition matrix of the 
filter. Hence, it is possible that this added knowledge will allow the filter to estimate a 
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system with more unknowns than measurements. Based on the developments in the 
section above, the carrier phase observation matrix for the seven satellite case is 
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where the parenthetical superscript refers to the GPS satellite, h is the unit vector to the 
SV, (jt,y,z) is the receiver position, ( x,y,z ) is the receiver velocity, and N U ’ J) is the 
integer ambiguity difference between the i th and j th SVs. The methods for utilizing this in 
propagating the Kalman filter are the same as in Section 2.7, The Kalman Filter. 


3.3.1. Starting Requirements 

At the end of the homing phase, the position, velocity, and the associated 
covariances remain the same; the integer ambiguity and an initial value for its covariance 
must be established. The true initial integer ambiguity difference is calculated from the 
true state (r) and the positions of the two SVs in question (p,, p ( ) 


N u ' n = fix 


R“ F l|-||p7- f ||) /x 


(56) 


where fix() is the function that rounds the argument towards zero, and X is the wavelength 
of the GPS carrier. The initial estimate of the integer ambiguity is calculated from this 
same formula, but replacing the state estimate (r ) for the true state. 

The integer ambiguity covariance matrix is calculated from the position 
knowledge at the end of the homing phase. If the position is known from the homing 
phase to a x , a y , and a z , then this establishes the appropriate range for the integer 
ambiguity uncertainty as well, due to their relation through the carrier wavelength. The 
trace of the covariance matrix is used because this will yield a conservative value. Also, 
since the integers to be estimated are actually integer differences, the range of the 
uncertainty is increased by a factor of ~J2 . Thus, the integer ambiguity standard deviation 
is given by 
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(57) 


a N =^V 2 ‘ Trace(I V 

The state noise for the positions and velocities were computed the same way as in 
the homing phase. As described in Section 2.7.2, 10% of the acceleration due to drag at 
the rendezvous altitude was used. The integer ambiguity state noise was determined from 
the position state noise using the factor -J2 / X, as above. 

The measurement noise for the phase double differences is determined first by 
establishing the knowledge on a single phase measurement. This study investigates two 
different knowledge levels on this measurement, 15° and 5°; this will provide an idea of 
the measurement accuracy required. The 15° level is the worst case phase knowledge 
measured in a ground test of single difference GPS attitude determination 29 , and 5° is 
used as a more optimistic estimate. Because the actual measurement used in the filter is a 
double difference, a factor of V2 ■ V2 = 2 is utilized to account for the differencing 
process. For example, if it is assumed that the measurement noise level is 15°, then the 
double difference noise will be on the order of 30°. 

3.3.2. Correlated Phase Measurement Errors 

As discussed in Section 2.7.3, Correlated Pseudorange Measurement Errors, an 
assumption in the Kalman filter is that the measurement noise is uncorrelated. As with 
the pseudorange, the phase sampling rate cannot be so high that the measurement noise 
becomes correlated. A phase sampling interval of two seconds has been suggested in 
prior discussions 30 , so this is the value used in the current investigation. 

3.4. GPS Satellite Tracking 

At the start of the .stationkeeping phase, it is necessary to select seven satellites for 
use in the resolution of the integer ambiguities. This is accomplished by propagating the 
satellite orbits and the receiver orbits, and determining line of sight for all the SVs until 
they drop out of view. In this simulation, this is an easy task because the orbits are 
already calculated and stored in tables. In a real mission, this is a solvable problem as 
well, because all satellite orbits are predictable. The seven satellites are then selected 
simply by choosing those that will stay in view the longest. 

There is the question of whether or not these seven satellites will stay in view for 
the entire process of resolving the integer ambiguities. In this simulation, the resolution 
process took on average about 10 minutes, and it was possible to maintain seven satellites 
in view for this period of time. If, however, a satellite were to drop out of view before the 
integers were resolved, the solution would be to start the integer resolution process over 
again with seven new satellites and resolve the new ambiguities. With the optimal 21 
satellite constellation in place, it will not take long for an acceptable set of satellites to 
come into view. This solution is plausible for all the missions this investigation is 
concerned with, because the delay of a few orbits will have no significant impact on 
accomplishing the mission objectives. As mentioned previously, the velocity uncertainty 
in the chaser (2- 10 2 m/s) is not critical; even if the resolution process takes an hour, this 


41 



velocity maps into a displacement of less than 100 m. The chaser will still be far enough 
from the target to ensure safety, while not imposing restrictive AV requirements. 

3.5. Stationkeeping Results 

3.5.1. State Knowledge: 15° Measurement Noise 

Figures 3.5. l.a-c show the knowledge of the position, velocity, and ambiguities 
for a standard stationkeeping phase, assuming 15° measurement noise. In each case, the 
advantages of differential GPS are shown; the knowledge benefits are two orders of 
magnitude or more over standard C/A positioning. The position knowledge decreases 
from around 20 m to 0.03 m, the velocity from 2- 10' 2 m/s to 2- 10' 5 m/s, and the 
ambiguities from 30 wavelengths to 0.2 wavelengths. 



Figure 3.5.1. a: Stationkeeping Position Knowledge and Residual, 
15° Measurement Noise 



Figure 3.5. l.b: Stationkeeping Velocity Knowledge and Residual, 
1 5° Measurement Noise 
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Simulation Time (min) 

Figure 3.5. l.c: Stationkeeping Ambiguity Knowledge and Residual, 

15° Measurement Noise 

An interesting point to note is that while a static survey on the ground (no receiver 
motion) typically requires 20-30 minutes to completely resolve the integer ambiguities, 
the kinematic method for the satellites in orbit only requires about 10 minutes. This may 
be due to the observability of the systems in question. For the static initialization, the 
resolution of the ambiguities is almost completely dependent on the motion of the GPS 
satellites to create a sufficiently observable system. However, for the case of the satellites 
in orbit, the position of the receiver is changing much more rapidly, possibly resulting in 
increased observability of the SVs, and hence a faster convergence time. 

Another result to notice is that the ambiguities are resolved in the orbital case in 
less than half the time than a ground use of Kinematic on the Fly 31 . Most likely, this is 
due to the fact that the mechanics of the system are exploited in the rendezvous. If the 
state is simply allowed to random walk, then there is no knowledge of the relation 
between the states at one epoch and the next. However, the mechanics in the rendezvous 
are well known and thus can be used to enhance the knowledge of the state. The overall 
effect, then, is to provide faster solution convergence. 

3.5.2. State Knowledge: 5° Measurement Noise 

The results of using more accurate measurements is shown in Figures 3.5.2.a-c. 
One result is clear: the state knowledge converges at a faster rate. Whereas the above 
poorer measurement case required approximately 10 minutes to converge, this case only 
requires about 7 minutes. This observation is purely academic. Because the chaser can 
generally remain in stationkeeping mode until the ambiguities are known, the time to 
convergence is not an important factor here. 


43 



V 



Figure 3. 5. 2. a: Stationkeeping Position Knowledge and Residual, 
5° Measurement Noise 



Figure 3.5.2.b: Stationkeeping Velocity Knowledge and Residual, 
5° Measurement Noise 



Figure 3.5.2.c: Stationkeeping Ambiguity Knowledge and Residual, 
5° Measurement Noise 
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3.5.3. GPS Satellite Usage: 3, 5, and 7 SVs 

The results in the above sections were obtained through the use of 7 GPS 
satellites. However, as described in Section 3.3, Kinematic on the Fly GPS and Kalman 
Filtering , it may not be necessary to utilize this full set. Because the system dynamics are 
known, this added knowledge may make it possible to achieve comparable accuracies 
while using fewer SVs. Note that by tracking fewer satellites, fewer ambiguities need to 
be estimated in the Kalman filter. For this part of the investigation, it is assumed that the 
measurement noise is 5°. Figures 3.5.3.a-c show that it is indeed possible to utilize only 
5 satellites while still achieving results on the same order as the 7 satellite case. 



Simulation Time (min) 

Figure 3. 5. 3. a: Stationkeeping Position Knowledge and Residual, 5 SVs 



Figure 3.5.3.b: Stationkeeping Velocity Knowledge and Residual, 5 SVs 
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Figure 3.5.3.c: Stationkeeping Ambiguity Knowledge and Residual, 5 SVs 

When only 3 GPS satellites are used in the state estimation (again with 5° 
measurement noise), the filter performance is clearly degraded. The convergence time is 
more than twice as long as when 5 or 7 SVs are used, and even then. Figures 3.5.3.d-f 
show that the entire state knowledge (position, velocity, and ambiguity) is an order of 
magnitude worse than the other cases. Thus, it is not recommended to utilize fewer than 
5 GPS SVs in the differential positioning, unless 1 m uncertainty is acceptable. 



Figure 3.5.3.d: Stationkeeping Position Knowledge and Residual, 3 SVs 
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Figure 3.5.3.e: Stationkeeping Velocity Knowledge and Residual, 3 SVs 



Figure 3.5.3.f: Stationkeeping Ambiguity Knowledge and Residual, 3 SVs 
3.6. Stationkeeping Conclusions 

This investigation of the stationkeeping phase of the rendezvous has led to several 
important results. First, Kinematic on the Fly GPS can be used to estimate a spacecraft's 
position, velocity, and integer ambiguities. Even though a twelve state system cannot be 
solved with seven satellites alone, the system dynamics provide additional information so 
that the procedure can converge. Second, the process of resolving the integer ambiguities 
takes about 10 minutes, less than half the time of a ground resolution. Third, using as few 
as 5 SVs, the filter will converge, yielding position knowledge on the order of a 
decimeter. Fourth, a decrease of the level of measurement noise from 15° to 5° translates 
into to about a three minute reduction in convergence time. Finally, obtainable final 
position, velocity, and ambiguity knowledge values are 0.03 m, 2- 10 5 m/s, and 0.2 
wavelengths, respectively. 
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4. Terminal Approach via Differential GPS 

4.1. Overview 

The analysis of the stationkeeping phase showed that by maintaining a quasi- 
constant distance, the chaser and target relative positions could be determined to within a 
decimeter. Best results are obtained when utilizing a 7 satellite model and 5° 
measurement noise. For the terminal investigation, the 7 satellite method is used also. 
However, a conservative value of 15° measurement noise is used, partially to account for 
multipath errors; this is in accordance with several previous studies 32 . Once the integer 
difference ambiguities are resolved, the positions are known, and the chaser may perform 
a maneuver to rendezvous with the target. As with all burns, this AV is perturbed to 
account for errors in magnitude and direction. This places the chaser on a slightly 
incorrect trajectory, so that one or more terminal corrective burns will be required. 

The Kinematic on the Fly GPS method is utilized as in the stationkeeping phase. 
Given that the results from the stationkeeping phase prove that the "Kinematic on the Fly" 
method works for the orbital system, the main question to be answered in this study is 
how robust this method is to cycle slips. Recall that this occurs when the receiver 
miscounts the carrier phase, which in turn leads to errors in the state estimate. This was 
not a critical topic in the stationkeeping analysis, because the target and chaser are not on 
an approach trajectory. Thus, if a cycle slip occurred, the ambiguity resolution procedure 
could be restarted without risking mission success. However, when the chaser has 
entered the terminal phase, it is necessary to have access to the higher accuracy of carrier 
phase information. If a cycle slip occurs, there may not be enough time remaining to 
resolve the new ambiguities. In this case, the only information available would be from 
C/A pseudorange measurements, but position uncertainty on the order of C/A 
pseudorange errors is unacceptable for rendezvous. One possible solution would be to 
perform a burn to stop all relative motion between the chaser and target, and to enter 
another stationkeeping mode to resolve the ambiguities. However, it is desired to 
investigate whether cycle slips can be repaired in real time; these issues, as well as typical 
simulation results, are discussed in the following sections. 

4.2. GPS Satellite Tracking 

The terminal phase is assumed to last for 45 minutes, so a constellation change 
will have to be made at some point in this portion of the rendezvous. The method for 
selecting a different set of SVs is the same as in the stationkeeping phase, described in 
Section 3.4, GPS Satellite Tracking. In addition, the same method of estimating the 
ambiguity and the corresponding covariance is used as in the previous phase, outlined in 
Section 3.3.1, Starting Requirements. It is important to realize that because the 
ambiguities were resolved in the stationkeeping phase, the chaser position is known to 
less than one wavelength. Therefore, when a new set of SVs is chosen, the ambiguities 
are basically known immediately. This will allow the chaser to continue along the 
trajectory to the target, with complete state knowledge, no matter how many times the 
selected satellites are changed. 
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4.3. Kinematic on the Fly GPS and Terminal Filtering 

The same filtering procedure is utilized in the terminal phase as in the 
stationkeeping phase, described in Section 3.3, Kinematic on the Fly GPS and Kalman 
Filtering. In the stationkeeping phase, the main goal of the filtering was to resolve the 
integer ambiguities; this yields the relative positions needed to begin the terminal phase 
processing. Once the integers are known, one might argue that it is only necessary to 
filter the position and velocity for the terminal phase, because this represents all the 
needed information. However, it may be possible to utilize the Kalman filter to detect 
and fix cycle slips; thus, the integer states are filtered as well. Even if this method cannot 
be used to repair the slips, their detection alone could potentially save the mission. 
Knowledge that a cycle slip has occurred permits the chaser to perform a bum to enter 
another stationkeeping mode. In this way, the state knowledge can be reestablished, and 
the rendezvous continued. 

The same procedures to determine the filter starting requirements are used in the 
terminal phase as in the stationkeeping phase. The position and velocity knowledge and 
states remain the same. The integer ambiguity and its covariance are determined through 
the procedure described in Section 3.4, GPS Satellite Tracking. The state and 
measurement noise levels are identical to those in the stationkeeping phase. 

4.4. Corrective Maneuvers 

The procedure for the midcourse burns in the homing phase is utilized here as 
well. The AV is assumed to be instantaneous, so that only the velocity state and 
covariance change. The velocity estimate after the burn is set to the nominal velocity, and 
the true velocity takes into account perturbations in the magnitude and direction of the 
burn. Upon investigating several simulations, it turns out that the terminal burn errors are 
extremely small, for two reasons. First, the chaser trails the target by only 2 km, so the 
initial burn is very small (6 10 4 m/s), and in turn, the initial burn errors are very small. In 
addition, the position is known so well at the end of the stationkeeping phase, that errors 
due to position uncertainty are small. Typical midcourse corrections at about 10 minutes 
before rendezvous are only on the order of 10 4 m/s. As a result, it was not deemed 
necessary to attempt to develop a comprehensive guidance strategy for the terminal 
phase; the minuscule propellant savings would not be worth it. A single maneuver about 
10 minutes before rendezvous should suffice. 

4.5. Terminal Phase Results 

4.5.1. State Knowledge, No Cycle Slips 

This section highlights the state knowledge obtainable in the terminal phase of the 
rendezvous. The results improve only slightly over the stationkeeping phase, most likely 
because of the increase in the amount of processed data. Figures 4.5. l.a-c show the state 
knowledge for a typical terminal phase. The position knowledge improves to around 1.3 
cm, while the velocity knowledge remains around 10 5 m/s. The integer ambiguity 
knowledge increases to around 6- 10 2 wavelengths. The times of constellation change are 
seen most clearly in Figure 4.5. 1 .c, where they are shown by small discontinuities in the 
dotted (covariance) line at approximately 123, 130, and 139 minutes. The most 
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important aspect of these results is that the state knowledge can be maintained even 
though constellations are changed. 



Simulation Time (min) 

Figure 4.5.1. a: Terminal Position Knowledge and Residual 



Simulation Time (min) 

Figure 4.5. 1 .b: Terminal Velocity Knowledge and Residual 
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Figure 4.5. l.c: Terminal Integer Ambiguity Knowledge and Residual 
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4.5.2. Cycle Slips 

A potentially serious event when using carrier phase relative positioning is the 
occurrence of a cycle slip. These are caused when the receiver zero crossing counter 
either misses one or more crossings, or mistakenly increments the counter. A cycle 
increment could be missed if there is some type of drop in signal strength; whereas an 
electromagnetic surge may cause the counter to increment too many cycles. The end 
result is that instead of position knowledge on the order of centimeters, the error is now at 
least one wavelength (19 cm), and possibly several wavelengths. The effect of such an 
occurrence on position knowledge is shown in Figure 4.5.2.a, where a cycle slip of one 
wavelength has been artificially inserted. Of course, it is possible that many cycles will 
be skipped, and as the number of skipped cycles increase, so does the position error. 



Simulation Time (min) 

Figure 4. 5. 2. a: Terminal Position Knowledge and Residual, Cycle Slip at 1 18 Minutes 

Unfortunately, the state knowledge has become so small at the time of the cycle 
slip that the Kalman filter basically rejects the information contained in the measurement. 
As a result, the filter cannot process such an occurrence on the fly, and additional logic is 
required. The key to detecting cycle slips lies in the examination of the measurement 
residual, given by 

Az = z - Hx (58) 

The statistics of the measurement residual should be consistent with the measurement 
noise levels used in the Kalman filter. Consider, for example, the measurement residual 
time history for the case above, shown in Figure 4.5.2.b. 
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1.2 



Simulation Time (min) 


Figure 4.5.2.b: Terminal Phase Measurement Residual 

The measurement residuals up to the time of the cycle slip appear to be stationary, 
and they are consistent with the Kalman filter measurement noise levels. Recall that the 
measurement noise was set at 30°, or 0.0833 wavelengths; this is the standard deviation 
of the measurement residuals in Figure 4.5.2.b. At the time of the cycle slip, the residual 
makes a significant jump, and is not at all consistent with the filter statistics. In the 
simulation program, the flag for a cycle slip is whenever a measurement residual lies 
outside the 4o value. The choice of this value allows the simulation to ignore consistent 
measurement residuals outside the 3c r value, while flagging statistically inconsistent 
points at times of cycle slips. 

An advantage of this detection method is that the measurement residual basically 
jumps to the number of cycles that have been skipped. Note that the residual does not 
jump exactly to the skipped number, because state and measurement noise is added to the 
system. However, simply rounding to the nearest integer should yield the correct number 
of skipped cycles, because the noise levels are two orders of magnitude less. Once this 
number is established, it is a simple matter to repair the cycle slip. Recall that the integer 
ambiguities are actually integer differences, so the cycle slip could have occurred on the 
signal from either SV. Both possibilities are checked by taking the state before the cycle 
slip, and propagating the Kalman filter for the one time step. The case that yields the 
minimum measurement residual is chosen to be truth, and the processing is continued. 

The effectiveness of this procedure is highlighted in Figures 4.5.2.c-e. These 
show the position, integer ambiguity, and measurement residual time histories for the 
case above, with a cycle slip at 1 15 minutes. The repair scheme completely eliminates 
any effect of the cycle slip, thus permitting the completion of the mission. 
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Figure 4.5.2.c: Terminal Position Knowledge and Residual, Repaired Cycle Slip 
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Figure 4.5. 2.d: Terminal Integer Ambiguity Knowledge and Residual, 

Repaired Cycle Slip 
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Figure 4.5. 2.e: Terminal Measurement Residual, Repaired Cycle Slip 

While this method was demonstrated to only one integer state and using a slip of 
only one cycle, it can be extended easily to apply to all states, and it is effective for any 
number of skipped cycles. If additional logic is incorporated, this type of method would 
permit a rendezvous if several or all of the integer difference states skipped any number 
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of cycles at any number of times. Consider for example the case where both phase 
counters contributing to the same integer ambiguity difference experience a cycle slip 
during the same fdtering step. Because the measurements are phase differences, then at 
least one of the slips would show up in two separate measurements. From this second 
measurement, one of the cycle slips could be quantified, so that the other slip could then 
be determined from this new information and the first measurement. If these steps were 
added to the rendezvous program, then the mission would not be at risk during such an 
event. The possible scenarios involving the combinations of cycle slips in phase 
measurements become more and more complicated. However, the logic can be 
developed that would greatly reduce the chance of failure. 

4.6. Terminal Phase Conclusions 

Two major results can be reported as a result of the terminal phase investigation. 
First, the Kinematic on the Fly GPS method can be used for the terminal phase of the 
rendezvous, even though the selected SV constellation must be changed several times. 
Position and velocity knowledge on the order of 1 cm and 10 -5 m/s respectively can be 
achieved. In addition, a method of detecting and repairing cycle slips has been 
demonstrated. Only a few of the many possible cycle slip scenarios have been 
investigated; thus, more research is required before the repair method can actually be put 
to use. Regardless, the detection method is valid for all cases, and is an important result 
in itself. If a cycle slip occurs, and it is at least detected, then a maneuver can be 
performed to place the chaser in a stationkeeping mode, thereby reducing mission risk. 


5. GPS Attitude Determination 
5.1. Overview 

The investigation to this point has revealed that knowledge of the chaser position 
can be maintained at the centimeter level, making a rendezvous possible from this 
viewpoint. In addition to position knowledge, the vehicle attitude must be known 
accurately as well. Poor attitude knowledge can just as easily lead to mission failure as 
poor position knowledge. As a result, this investigation includes an analysis on utilizing 
GPS to determine vehicle attitude; however, an attitude control system has not been 
modeled. GPS attitude determination utilizes the carrier phase observable, similar to the 
position determination system. However, the integer ambiguity state is not needed, 
because of the close proximity of the receivers. The following "classic" assumptions are 
made in the procedure described below 33 

• the spacecraft is rigid 

• motions about the center of mass do not affect the motion of the center of mass 

• the body frame defines the principal axes 

• small rotations 

• near circular orbit 

• other external torques (aero, solar) are neglected 
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5.2. Attitude Equations of Motion 

The development of the attitude mechanics for a gravity gradient stabilized 
satellite basically follows that presented in several references 34 . An important difference 
is the reference frame convention; for this simulation, the nominal (no rotation) body 
frame is coincident with the instantaneous target orbital frame. Using a 3-2-1 Euler 
rotation (in this case, pitch - yaw - roll), the linearized equations of motion for the 
principal axes are 


I P P = ~ 3(/ r — I y )n 2 p + M p 
I y y + (7 r + I y - I p )nr + (/„ - I r )n 2 y = M y 


(59) 


i;r -(/,+/,- /„ )ny + (/, - /, )n 2 r = -3(1 p - I y )n 2 r + M r 


where I pyr are the, moments of inertia, n is the orbital angular velocity (mean motion), and 
M pyr are externally applied moments. For this investigation, the values for the moments 
of inertia are 


/ = 20 kg m 2 / = 6 kg • m 2 / r = 24 kg • m 2 


Neglecting external torques, the above equations in matrix form are 
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(61) 


The solution to this set of equations is obtained through the use of the matrix exponential 
method. Given a set of differential equations, 


x = Ax 


(62) 


the solution x(t) is found from 


x(f) = e A ' 


(63) 
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For the Kalman filter, the attitude state transition matrix is found using this method. 
Knowing that the measurements are taken at intervals of At, the state transition matrix in 
between measurements is simply 


0 = <? AA/ (64) 

Note that because A and At are fixed, the state transition matrix need be computed only 
once. 

5.3. Carrier Phase Measurements 

This brief discussion of carrier phase measurements is based on that found in 
Melvin and Hope 35 . Given a pair of GPS antennas, the phase measured at each is given 
by 


<t>. = k, a, - tor 

1 1 ' ( 6 ! 

<t> 2 = k 2 a 2 - cor 

where kj is the wave vector from the GPS satellite in question, a s is the position of the 
antenna phase center, to is the circular frequency of the GPS carrier, and t is the time of 
the measurement. This method carries with it two important assumptions; both antennas 
use the same oscillator to determine the phase, and the measurements are made at the 
same time. In addition, if the plane wave approximation is made, and the antenna 
distance is assumed negligible compared to the distance to the GPS satellite, then the 
wave vector is given by 


k, = k = -K p = — p (66) 

A 

where k=2k/X is the wave number, X is the carrier wavelength (approximately 19.0 cm 
for the LI frequency), and p is the unit vector from the receiver towards the GPS 
satellite. A phase difference is constructed from these two measurements 

A<|> = <J> 2 — <j)j = k (a 2 -a,) = k b (67) 

where b = (a 2 — a,) is the baseline vector between the antennas. These phase differences 
are the measurements from which attitude will be determined, as described in the section 
below. 

5.4. Attitude Determination from Carrier Phase Measurements 

In order to determine vehicle attitude, this development assumes that there is a 
known nominal attitude that the spacecraft would traverse in the absence of perturbing 
forces. In the case of gravity gradient stabilization, the chaser is nominally aligned with 
the target orbital frame. The true attitude is then this nominal attitude plus a small 
deviation resulting from perturbing forces and errors. The simulation utilizes this 
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deviation in pitch, yaw, and roll (and their time derivatives) as the attitude state. 
Similarly, the measurements are deviations in the phase differences from the true and 
expected states. In effect, the double difference measurement is used in the attitude as 
well as the position determination. The only difference is that in the position 
determination, the double differences were between receivers and then satellites , whereas 
in the attitude determination, the differences are between the receivers and then the true 
and estimated state. The model described below is based on the developments in several 
references 36 - 37 . 

Considering the same pair of GPS antennas above, with coordinates from the 
center of mass a, and a 2 expressed in the body frame, the baseline vector in the body 
frame, b b , is 


b b = a 2 -a, (68) 

This baseline is assumed to be known to less than A/100 from standard measurement 
techniques, including possibly static GPS measurements. This vector can be transformed 
to the target orbital frame via the direction cosine matrix cp 

b [of = cpb b (69) 

The direction cosines are assumed to be composed of a nominal and a perturbed 
component, so that 


<P = Sep cp 0 (70) 

where cp 0 is the direction cosine matrix for the nominal attitude, and 8(p is the direction 
cosine matrix using the small angle approximation. For the 3-2-1 Euler sequence, these 
matrices are 


and 


<Po 


CyCp CySp —Sy 

-CrSp + SrSyCp CrCp + SrSySp SK2y 
SrSp + CrSyCp —SrCp + CrSySp CrQy 


Sep = 


1 8 p -5 y 

-8 p 1 8 r 

8y — 8r 1 


( 71 ) 


(72) 


where S and C are sine and cosine, respectively. Upon evaluating <p = Sep cp 0 , the result 
can be written in the form 


57 



tp = tPo + A(p 


(73) 


where 


A<p n = (SrSyCp - CrSp)bp - ( SrSp + CrSyCp)by 

A<p )2 = ( CrCp + StiySp)bp - (CrSyS/7 - SrCp)by 

Acp 13 = SrCybp - CrCyby 

Acp 2 i = —CyCpbp + (SrSp + CrSyCp)br 

Acp 22 = -CySpbp + (CrSySp — SrCp)br 

A(p 23 = Sybp + CrCybr 

Acp 3l = CyCpby -(SrSyCp-CrSp)br 

Acp 32 = CySpby — (CrCp + SrSySp)br 

Acp 33 = -Sy5y - SrCybr 


Introducing the small angle approximation into the A(p matrix yields 


Acp = 
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-5 P 
5y 


bp -by 
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(74) 


(75) 


The baseline measurement in the target orbital frame can be found from 

b tof = (<Po + A( P) b b (76) 

This result can be applied to the phase difference equation to yield 

2.K — _ 2i t — 2k * . ■ 

A<J) = k b tof = — — p (cp 0 + A(p)b b =-Y P <Po-y P' Atpb b (77) 

As stated above, a nominal attitude state is assumed; corresponding to this, there will be a 
nominal phase difference, A<b 0 , and a differential phase difference, 5A<{>, so that the total 
phase difference can be expressed as 

A<j> = A({) 0 + 5A(J) (78) 


where 


, i 2 k _ 

A(t>0= "T p (p0 


(79) 


58 



5A* = - t P-A<pb b 


Because the nominal attitude is known for the entire trajectory, the nominal phase 
difference is known as well; hence, the quantity of interest reduces to the differential 
phase difference 5A(J). 

In order to determine how the measurement 5A<)> relates to the state variables 
P,y,r,p,y,f, it is necessary to determine the observation matrix. This is done by simply 
evaluating the above equation 
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Finally, the observation matrix is 
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or, more succinctly. 
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This equation represents a single phase measurement. The incorporation of more 
measurements is done simply by appending rows to the observation matrix, so the final 
result looks something like 
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This is the measurement equation that is used in the attitude Kalman filter, outlined 
below. 

It is important to remember that a critical assumption in this formulation is that 
the attitude deviation from nominal is small (less than A/b) 38 . So long as this criterion is 
met, then the differential phase measurement completely defines the state of the system, 
and it is not necessary to take into account the integer ambiguity. However, if this 
assumption is violated, then the integer ambiguity must be resolved before the attitude 
can be positively determined. While this is certainly possible, it would increase the 
complexity of the model. 

5.5. The Attitude Kalman Filter 

The attitude Kalman filter utilizes the same procedure outlined in the previous 
developments. The attitude state is the deviation from nominal of the spacecraft pitch, 
yaw, roll, and their time derivatives 

x = [p y r p y rf (84) 


The state transition matrix is determined from the Euler 3-2-1 linearized equations of 
motion, as described in Section 5.2, Attitude Equations of Motion. Also, the observation 
matrix for the measurement equation is derived in Section 5 A, Attitude Determination 
from Carrier Phase Measurements. Otherwise, the filter is propagated utilizing the 
techniques described in Section 2.7.1, Filtering Procedure. 

Several parameters must be specified to initialize the filter. A standard 
assumption is that all covariance and noise matrices are initially diagonal. The initial 
attitude knowledge of the chaser (the covariance matrix P 0 ) must first be established. 
Assuming only coarse attitude control using gravity gradient until the rendezvous is 
underway, a knowledge level of 10° is assumed 39 . The nominal state is zero attitude, so 
the initial truth and estimate are perturbed around this state consistent with the given 
initial covariance. 

The main contribution to the state noise is assumed to be torques due to errors in 
moments of inertia. The torque due to gravity gradient effects is 40 



( 85 ) 


Combining this with the relation 
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( 86 ) 


d 2 e _ T 
dt 2 ~ / 


yields appropriate values for the standard deviations of the state noise and state rate noise 
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where / 9 and / Q are scale factors, typically set to 0. 1 to account for errors in system 
modeling. 

Because the measurement for the attitude filter is exactly the same type as for the 
differential filter, the same measurement noise matrix is used. As described in Section 
3.3.1, Starting Requirements, this consists of a diagonal matrix consisting of either 15° or 
5°, which is then corrected for the double differencing process by a factor of 2. 


5.6. Attitude Results 

5.6.1. Two Receiver Investigation 

Given two GPS receivers and no other information, it is possible to determine 
only two of the angles describing the baseline attitude. The third angle, the rotation about 
the baseline vector, cannot be established. However, once the information inherent in the 
attitude equations of motion is applied, it is possible to estimate the third unknown 
angle 41 . For the linearized attitude equations, the roll and yaw axes are coupled; thus, 
measurements corresponding to one axis will be reflected in the other through the system 
dynamics. As a result, a baseline aligned along either the roll or yaw axis will provide 
information on all three angles. Consider Figures 5.6. l.a-f, which illustrate the attitude 
knowledge assuming a phase measurement capacity of 15°, and measurements from 4 
GPS SVs. For each case, the knowledge of the observable axes starts around 1° for the 
homing phase, and decreases to about 0.1° for the terminal phase. For the unobservable 
axis, the knowledge seems to depend on which baseline was used. For the roll aligned 
antenna, the knowledge starts around 200° and falls to about 10°; the yaw aligned antenna 
begins at the same level, but falls to around 1°. For both cases, this poor knowledge 
eliminates the use of only two receivers for an attitude control system. 
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Figure 5.6.1. a: Homing Attitude Knowledge and Residual, Roll Aligned Antenna 



Figure 5.6. l.b: Stationkeeping Attitude Knowledge and Residual, Roll Aligned Antenna 



Figure 5.6. l.c: Terminal Attitude Knowledge and Residual, Roll Aligned Antenna 
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Simulation Time (min) 

Figure 5.6. l.d: Homing Attitude Knowledge and Residual, Yaw Aligned Antenna 



Simulation Time (min) 

Figure 5.6. l.e: Stationkeeping Attitude Knowledge and Residual, Yaw Aligned Antenna 



Simulation Time (min) 

Figure 5.6. l.f: Terminal Attitude Knowledge and Residual, Yaw Aligned Antenna 

Unlike the other two axes, the pitch angle is not coupled with any other motion. 
Hence, if the baseline vector is oriented along the pitch axis, it is impossible to observe 
motion about that axis. This is clear in Figures 5.6.1.g-i; the pitch axis is completely 
inestimable. However, the knowledge of the remaining axes falls from about 1° in the 
homing phase to about 0. 1° in the terminal phase. An interesting point to note is the 
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oscillatory nature of the pitch residual; this phenomenon is consistent with the pitch 
period. 



Figure 5.6. l.g: Homing Attitude Knowledge and Residual, Pitch Aligned Antenna 



Figure 5.6. 1 .h: Stationkeeping Attitude Knowledge and Residual, Pitch Aligned Antenna 



Figure 5.6. l.i: Terminal Attitude Knowledge and Residual, Pitch Aligned Antenna 
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5. 6. 1.1. The Use of More GPS Satellites 

All of the above results utilize measurements from four GPS satellites. It is 
conceivable that more measurements could improve the attitude knowledge, so the effect 
of utilizing 6 SVs was investigated. Figure 5. 6. 1.1. a shows that while there is a small 
gain in the observable directions, the unobservable axis is still too poor to utilize for a 
control system. While this figure only shows the stationkeeping phase of the rendezvous, 
the results are similar for the homing and terminal phases. 



Figure 5.6. 1.1. a: Stationkeeping Attitude Knowledge and Residual, 6 GPS SVs 

5.6. 1.2. Increasing the Phase Measurement Accuracy 

If the measurement noise levels are reduced from 15° to 5°, while using 4 SVs, 
the state estimates should of course increase in accuracy; the question remains whether 
the use of only a pair of GPS antennas will supply enough accuracy in all dimensions for 
use in a control system. Figure 5.6. 1 ,2.a-c show the results of this investigation. During 
the homing phase, the attitude knowledge in the unobservable direction still remains on 
the order of 10°; however, the knowledge of the remaining axes is established at about 
0.5°. In the stationkeeping and terminal phases, the accuracy of the observable axes falls 
to slightly below 0.1°. The largest increase is evident in the unobservable axis; its 
knowledge falls to the order of 1°. 
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Figure 5.6. 1.2. a: Homing Attitude Knowledge and Residual, 5° Measurement Noise 



Simulation Time (min) 


Figure 5.6. 1 ,2.b: Stationkeeping Attitude Knowledge and Residual, 
5° Measurement Noise 



Figure 5.6.1.2.c: Terminal Attitude Knowledge and Residual, 5° Measurement Noise 
5.6.2. Three Receiver Investigation 

In this section, the use of three GPS antennas is considered. In this way, more 
than one baseline vector is available, thus eliminating the problem of the unobservable 
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angle seen in the above sections. Three GPS SVs are utilized, providing three double 
differences across two different baselines for a total of six measurements. 

5.6.2. 1. Phase Measurement Capacity: 15° 

For even the worst case measurement noise values, the results are very promising, 
as shown in Figures 5.6.2. l.a-c. The knowledge of all three attitude angles is maintained 
on the order of 1° for the entire homing phase. The estimates improve in the 
stationkeeping and terminal phases, to around 0.1°, three axis. 



Figure 5.6.2. 1. a: Homing Attitude Knowledge and Residual, Three Antenna 



Simulation Time (min) 

Figure 5.6.2. 1 .b: Stationkeeping Attitude Knowledge and Residual, Three Antenna 
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Figure 5.6.2. l.c: Terminal Attitude Knowledge and Residual, Three Antenna 

5.6.2.2. Phase Measurement Capacity: 5° 

The decreased measurement noise directly affects the knowledge of the attitude 
state. As with the 15° case, all three states are estimated consistently due to the inclusion 
of the third antenna; this is shown in Figures 5.6.2.2.a-c. However, the lower levels of 
noise lead to homing phase attitude knowledge on the order of 0.3°, falling to about 0.05° 
in the terminal phase. 



Figure 5.6.2.2.a: Homing Attitude Knowledge and Residual, Three Antenna 
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Figure 5. 6.2. 2. b: Stationkeeping Attitude Knowledge and Residual, Three Antenna 



Figure 5.6.2.2.c: Terminal Attitude Knowledge and Residual, Three Antenna 

5.7. Attitude Determination Conclusions 

The main consequence of this investigation is that GPS can be utilized to 
determine the attitude of a spacecraft to within about 0.05°. Several other results can be 
reported as well. For example, it is possible to use only two GPS receivers, coupled with 
the attitude dynamics, to determine spacecraft attitude. The best case scenario was when 
the baseline was along the yaw axis, yielding knowledge on the order of 1°. However, 
this two receiver method is not recommended for use in an attitude control system, even 
if more SVs are used or the measurement noise level is decreased. Also, if three GPS 
receivers are utilized, then full attitude knowledge can be obtained. If a 15° measurement 
noise level is used, the attitude knowledge converges to about 0.1°; whereas, if a 5° 
measurement noise level is used, 0.05° is the result. 


6. Summary and Recommendations 

This study has investigated a wide variety of areas in the use of GPS for 
automated rendezvous and docking. The main result is that it is possible to utilize GPS 
alone for the position and attitude determination of such a mission. It is important to 
realize that a simulation of the actual docking procedure has not been created; thus it is 
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not appropriate to directly compare rendezvous control requirements and the knowledge 
achievable from GPS. However, the GPS state knowledge must be better than the control 
requirements, or else the GPS system cannot be used. Thus, assuming that the 
requirements for a rendezvous are those shown in Table 6. a, state knowledge from GPS 
has been shown to exceed the necessary values, consistent with the assumptions in this 
investigation. 

Table 6. a: Rendezvous Control Requirements and GPS Knowledge 


Lateral Misalignment 
Docking Speed 
Angular Offsets 

Of course, it is advisable to include redundant systems, perhaps other than GPS based, to 
ensure mission safety. This section will briefly highlight some other results, and outline 
subjects that require further research. 

This work has shown that the use of the Kalman filter is appropriate for the 
homing phase of a rendezvous, even in the presence of selective availability. State 
knowledge on the order of 20-30 m is obtainable using C/A code positioning. The 
Kalman filter was also demonstrated to be applicable to the resolution of the integer 
ambiguities in the stationkeeping mode, and to the accurate estimation of the state in the 
terminal phase. Position and velocity knowledge on the order of 1 cm and 10 5 m/s 
respectively can be achieved at rendezvous. In the homing phase, guidance strategies and 
the effects of mismodeling SA statistics were investigated. "Kinematic on the Fly" GPS 
was analyzed in the stationkeeping phase; it was proven to be applicable, even when 
using as few as 5 SVs. Cycle slips were examined in the terminal phase, and a method of 
at least detecting such a slip, if not repairing it, was demonstrated. Attitude knowledge 
on the order of 0.05° was shown to be achievable. 

While this research has laid a strong foundation for the use of GPS in a 
rendezvous mission, there still remain many areas that need to be investigated further. 

For example, this study assumed that the SA noise was stationary with zero mean; 
however, this may not always be true. Also, the use of a sequential filter may not be 
entirely appropriate when attempting to estimate integers; as mentioned previously, the 
use of a batch filter may be more advantageous. In the discussion of cycle slips, a method 
of detecting such an occurrence was shown. This knowledge allowed the cycle slip to be 
repaired, and the state knowledge to be maintained. However, a more in depth 
investigation into the logic to repair concurrent cycle slips (on one or more channels) 
must be performed. In addition, this research utilized larger phase measurement errors to 
compensate for possible multipath errors; if a more exact multipath model could be 
developed and utilized, system errors could be reduced. Further, the use of an "all in 
view" method of positioning could result in decreased error levels. Finally, the effects of 
uncertainties in the positions of the GPS SVs, perhaps due to a ground station failure, 
could be investigated. All of these areas and more require further research before the 
methods explored in this study can be implemented in an actual rendezvous mission. 


Control 

Requirements 1 42 
7.5 cm 
20 mm/s 
1.5° 


GPS Knowledge 
1 cm 

0.01 mm/s 
0.05° 
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Appendix A: Description of Computer Programs 


A. 1. Overview 

The computer programs used in this study were developed on an IBM 486/50DX 
machine. The code was written in the MATLAB 4.0 for Windows programming 
environment. Three main driver programs have been developed which address the needs 
of this investigation: 


1) SA - generates the selective availability data via the Braasch model 

2) RVLOOKUP - creates the lookup tables with the GPS SV state information 

3) GOMONTE - simulates the rendezvous, using GPS position and attitude 

determination, from homing through terminal phases 


These drivers call the necessary subroutines to perform whichever task is desired. An 
outline of the control flow for each program is included below; each successive 
indentation represents a deeper level of subroutines. All programs and subroutines are 
discussed in greater detail in the following section. 

Program 1 
SA 

Program 2 
RVLOOKUP 

KEPLER 

FGC 

FCTORIAL 

FGS 

FCTORIAL 

MDLOS 

ANGBET 

CROSS 

D2R 

R2D 

ORB2XYZ 

D2R 

DCM313 

R2D 


Program 3 
GOMONTE 

ATTSTMC 

CROSS 

D2R 

DETVEL 
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GMPARMS 

RVMONTE 

ATTOBSM 

GPSMEAS 

GPSEQNS 

PHAS1SAV 

PHAS2SAV 

PHASEND 

RNDVECMT 

RV_SK 

ATTOBSM 

ATTSTMC 

GPSMEAS 

GPSEQNS 

PHASESVS 

R2D 

RNDVECMT 

RVKALM 

RNDVECMT 

RVLOC2E 

CROSS 

STATNOIS 

STMATRIX 

TRMOBSM 

RVCSTAT 

STMATRIX 

RVDELTAV 

CROSS 

DETVEL 

RVDOPSUB 

GDOP4 

AREATRI 

VOLTETRA 

RVINIT 

RNDVECMT 

RVDELTAV 

CROSS 

DETVEL 

RVKALM 

RNDVECMT 

RVLOC2E 

CROSS 

RVSVLOS 

MDLOS 

ANGBET 


72 



CROSS 

D2R 

R2D 

RVTERM 

ATTOBSM 

GPSMEAS 

GPSEQNS 

PHASESVS 

R2D 

RNDVECMT 

RVDELTAV 

CROSS 

DETVEL 

RVKALM 

RNDVECMT 

RVLOC2E 

CROSS 

STATNOIS 

STMATRIX 

TRMOBSM 

STATNOIS 

STMATRIX 

A. 2. Alphabetical Subroutine List and Description 

This section lists alphabetically all the subroutines used in the rendezvous 
simulation; a brief summary of each is included as well. All input and output variables 
are listed and described. Where necessary, the dimensions of matrices are explained in 
parenthesis following the variable description; for example, the parenthesis (x, y, z x SV,, 
SV 2 , SV 3 , ...) indicate that the three rows of the matrix are x, y, and z positions, while the 
columns refer to the satellites SV,, SV 2 , and so forth. Unless otherwise mentioned, all 
routines use as standard units kilograms, kilometers, and seconds. 


ANGBET - This function calculates the angle between two vectors, in radians. 

Inputs: 

x, y - two 3x 1 vectors 
Outputs: 

z - the angle between x and y, in radians 

AREATRI - This function calculates the area of the triangle formed from three endpoints. 
Inputs: 

a, b, c - the 3D endpoints of the triangle 
Outputs: 

x - the area of the triangle 
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ATTOBSM - This function calculates the attitude observation matrix, for use in the 
attitude Kalman filter. 

Inputs: 

ant_dist - distance between the receiver antennas, km 
ant_vect - unit vector from one antenna to the other 
wavelen - phase wavelength to use, km (standard LI = 19.0 cm) 
bvec2sat - matrix of unit vectors to the 4 GPS SVs (x, y, z x SV |t SV 2 , SV 3 , 
SV 4 ) 

Outputs: 

att_hcalc - the attitude observation matrix 

ATTSTMC - This function calculates the attitude state transition matrix. 

Inputs: 

momjn - moments of inertia, kg/km 2 
dt - time between measurements, sec 
mu - earth gravitational parameter, kmVs 2 
sc_rad - radius of the target (rendezvous radius), km 
Outputs: 

stm - attitude state transition matrix 

CROSS - This function returns the cross product of two 3D vectors. 

Inputs: 

x, y - two 3x1 column vectors 
Outputs: 
z - xxy 

D2R - This function converts degrees to radians. 

Inputs: 

x - angle in degrees 
Outputs: 

y - angle in radians 

DCM313 - This function returns the direction cosine matrix for a 3-1-3 rotation. 

Inputs: 

angl, ang2, ang3 - Euler angles in radians 
Outputs: 
y-DCM 

DETVEL - This function determines the velocity necessary to arrive at the origin in the 
desired time. 

Inputs: 

pos - position vector, km 
w - orbital angular velocity of target, rad/s 
tott - total travel time, sec 
Outputs: 
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vel - the velocity vector necessary to perform the maneuver, km/s 

FCTORIAL - This function calculates the factorial of an integer. 

Inputs: 

x - integer 

Outputs: 
y - x! 

FGC - This function is used in the universal variables approach to orbit propagation; it is 
only called by the KEPLER function. 

Inputs: 

z 

Outputs: 

C(z) 

FGS - This function is used in the universal variables approach to orbit propagation; it is 
only called by the KEPLER function. 

Inputs: 

z 

Outputs: 

S(z) 

GDOP4 - Given the unit vectors from the chaser to the four selected SVs, this function 
calculates the associated GDOP. 

Inputs: 

a, b, c, d - the unit vectors from the chaser to the four selected SVs 

Outputs: 

x - the associated GDOP 

GMPARMS - This routine specifies the desired simulation parameters. These parameters 
include everything, from the elevation angle to the phase measurement noise level to 
the number of Monte Carlo runs. 

GOMONTE - This is the main rendezvous simulation driver program. It calls all the 
necessary functions to determine both position and attitude in all phases of the 
rendezvous; this includes line of sight calculations, triangulations, etc. It can perform 
either a single simulation or a multi-run Monte Carlo simulation. The user inputs are 
defined in the routine GMPARMS. 

GPSEQNS - This function solves the GPS system of equations: 4 equations in 4 
unknowns, i.e. position (x, y, z) and clock offset. 

Inputs: 

xyztges - initial guess of position and clock offset, km and sec 
svpos - positions of the four SVs, km (x, y, z x SV,, SV 2 , SV 3 , SV 4 ) 
svrng - pseudoranges to the four SVs 
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rngtol - tolerance to iterate RANGE part to (i.e., 1 part in le-6) 

Outputs: 

xyzt - solution to the 4 nonlinear GPS equations, km and sec 
meascov - the inv(A)*inv(A)' matrix necessary to transform the covariance 
properly in the Kalman filter 

GPSMEAS - This function returns the actual position measurement (x, y, z) calculated 
from GPS; this includes the effects of SA. 

Inputs: 

svsoln - the positions of the 4 GPS satellites to be used in triangulation, km (x, 
y, z x SV,, SV 2 , SV 3 , SV 4 ) 

rtgt - the position of the origin of the local frame (i.e., the target position), km 
vtgt - the velocity of the origin of the local frame (i.e., the target position), km 
x - the position of the chaser in local coordinates, km 
sa - the SA noise values for the 4 satellites 
rngtol - the position tolerance to iterate to 
ijk21oc - DCM from inertial to local coordinates 
loc2ijk - DCM from local to inertial coordinates 
Outputs: 

gpsz - the actual GPS measurement (with SA), km 

meascov - the inv(A)*inv(A)' matrix necessary to transform the covariance 
matrix properly for the Kalman filter 

bvec2sat - the unit vectors from the receiver to the satellites, for use in the 
attitude observation matrix 

KEPLER - This subroutine is a standard orbit propagator; it uses the universal variables 
approach developed in Bate, Mueller, and White. 

Inputs: 

rO - initial position vector, km (3x1) 
vO - initial velocity vector, km/s (3x1) 
dt - time to propagate orbit, sec 
mu - gravitational parameter, kmVsec 2 
proptol - tolerance to iterate to find final r,v 
propiter - maximum number of iterations 
Outputs: 

r - position vector after time dt, km 
v - velocity vector after time dt, km/sec 

MDLOS - Given the receiver position and the SV position, this subroutine determines 
whether or not a line of sight exists. 

Inputs: 

sc_vec - the position of the observer, km 
sv_vec - the position of the satellite you're trying to see, km 
elevmask - the elevation angle from the earth tangent, deg 
alt - the altitude of the observer, km 


76 



Outputs: 

los - 0 - not in view, 1 - in view 

ORB2XYZ - Given the orbital elements of a body, as well as the gravitational parameter, 
this routine calculates the corresponding position and velocity vectors. This 
subroutine accepts several different possible combinations of orbital elements; 
however, a flag must be set to indicate which set is being used. 

Inputs: 

inputmode - which orbital element set is being used 
mu - gravitational parameter, kmVsec 2 
r - radius vector, km 
v - velocity vector, km/sec 
Outputs: 

orbelem - orbital elements output matrix 
flag - orbit descriptor 

PHAS1SAV - This program saves and then removes values at the end of the homing 
phase that are not needed in subsequent phases, thereby saving computer memory 
space. 

PHAS2SAV - This program saves and then removes values at the end of the 
stationkeeping phase that are not needed in subsequent phases, thereby saving 
computer memory space. 

PHASEND - This program saves and then removes values at the end of the terminal 
phase that are not needed in subsequent phases, thereby saving computer memory 
space. 

PHASESVS - This routine calculates which seven SVs stay in sight the longest from the 
current time; it is used in the carrier phase positioning in the stationkeeping and 
terminal phases. 

Inputs: 

svlos - the line of sight matrix; 0-not in sight, 1-in sight (time x SV,, SV 2 , 
SV„SV 4 ,...) 

Outputs: 

Iong_svs - the SVs that stay in view the longest 

svjimes - the times from current that the SVs will drop out of sight, sec 

R2D - This function converts radians to degrees. 

Inputs: 

x - angle in radians 
Outputs: 

y - angle in degrees 
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RNDVECMT - This function returns an error vector consistent with the statistics of the 
covariance matrix p. 

Inputs: 

p - covariance matrix 
Outputs: 

randvect - error vector consistent with p 

RV-SK - This is the main simulation routine for the stationkeeping section of the 
rendezvous. Based on the results from the homing phase, this routine calls the 
necessary programs to calculate the integer ambiguities, propagate the Kalman filter, 
and so forth. At the completion of this phase, the driver for the terminal section is 
called (RVTERM). 

RVCSTAT - This function calculates a variety of statistics on the C/A positioning. 
Inputs: 

tott - total time in the homing phase, sec 
tremain - time remaining in the homing phase, sec 
w - angular frequency of the target orbit, rad/sec 
p - current covariance matrix 
phi - DCM from inertial to b-plane 
x - current state truth 
xhat - current state estimate 
Outputs: 

stats - the diagonals of the covariance projected to the b-plane time of arrival 
covdiag - diagonal elements of covariance matrix 
xprojstate - the true state projected to b-plane arrival 
projstate - the estimated state projected to b-plane arrival 
sumcov - the harmonic mean (cube root of the product of the diagonal) of the 
position 3x3 covariance, projected to the b-plane time of arrival 
sumcov3 - not used 

xrerr - the truth radial error at b plane arrival, km 
rerr - the estimated error at b plane arrival, km 
dist - the current estimated distance to the target, km 

RVDELTAV - This function simulates a AV maneuver. Based on the desired velocity 
increment, the new true and estimated velocities are calculated, as well as the new 
velocity covariance. 

Inputs: 

tott - the total time in the homing phase, sec 

tremain - the time remaining in the homing phase, sec 

xold - true state before the maneuver 

xhatold - estimated state before the maneuver 

pold - covariance matrix before the maneuver 

w - the orbital angular velocity, rad/sec 

burnmagpct - the percent error in the magnitude of the burn 
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burnxpct - the percent error in the cross directions of the burn 
mode - 1 - the maneuver is to put the chaser on a rendezvous trajectory with 
the target 

2 - the maneuver is to stop relative motion between chaser and target 
(i.e., at the beginning of the stationkeeping phase) 

Outputs: 

x - true state after the maneuver 
xhat - estimated state after the maneuver 
p - covariance matrix after the maneuver 
normdv - the magnitude of the AV 

RVDOPSUB - Given the SVs in sight, this function determines the sub optimal 
combination of four (via the Noe, Myers, and Wu method) to use in the position 
determination. 

Inputs: 

svlos - the row vector of true or false (1 or 0) line of sight 

svpos - the row vector of all SV positions, km (SVj , SV, , SV, , SV 2 , 

sv 2y , ...) 

rtgt - the coordinates of the origin, km 
x - estimate of chaser position in inertial coordinates, km 
ijk21oc - DCM between inertial and local coordinates 
Outputs: 

svs - the four SVs to use in position determination 

svsoln - the row vector of the positions of the selected four SVs, km (SV, , 
sv, y , SV ]Z , SV 2x , SV 2y , ...) 
gdop - the GDOP value for the selected four SVs 

RVINIT - This subroutine basically performs a variety of initialization procedures. All 
the necessary variables for the computer simulation are created, and some of the 
default values are calculated. 

RVKALM - This function is a generic Kalman filter propagator. It accepts a mode 
variable that allows for special treatment of individual cases, such as SA 
incorporation or cycle slip repair. 

Inputs: 

mode - 2=no SA, 3=SA, 4=Gaussian SA, 22=cycle slip repair, 999=attitude 
stm - state transition matrix 
x - truth state 

xhatold - last state estimate 
pold - last covariance matrix 
q - state noise covariance matrix 
r - measurement noise covariance matrix 
h - observation matrix 
gpsz - actual GPS measurements 
Outputs: 
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xhat - new state estimate 
z - measurements 
p - new covariance matrix 


RVLOC2E - This function determines the direction cosine matrices from local to inertial 
frames and vice versa. 

Inputs: 

rtgt - the position of the origin of the local frame (i.e., the target position), km 
vtgt - the velocity of the origin of the local frame (i.e., the target position), km 
Outputs: 

ijk21oc - DCM from inertial to local 
loc2ijk - DCM from local to inertial 

RVLOOKUP - This program creates the lookup tables for the GPS satellite positions and 
velocities, as well as the target orbit, if desired. These lookup tables are stored and 
then used by GOMONTE to avoid doing those calculations repeatedly. 

Inputs: 

orbelem - orbital elements of the target (p, e, i, Q, CD, D) 

dt - time increment, sec 

endtime - total time to generate tables for, sec 

svelem - the elements of the GPS SVs (SV,, SV 2 , SV 3 , x p, e, i, £2, CD, u) 
elevang - the desired elevation angle, deg 
calctgt - flag to calculate the target states or not 
Outputs: 

rtgt - positions of the target, km (time x x, y, z) 

vtgt - velocities of the target, km/s (time x x, y, z) 

svpos - positions of the SVs, km (time x SV,,, SV ly , SV |Z , SV 2x , SV 2y , 

SV*,...) 

svvel - velocities of the SVs, km/s (time x SV lx , SV ]y , SV lz , SV 2x , SV 2y , ...) 
svlos - line of sight from target to SV, 1 = los, 0 = no los (time x SV lx , SV ly , 
sv, z , SV 2x , SV 2y , SV 2z , ...) 

RVMQNTE - This program performs one rendezvous simulation, from the homing 
through the terminal phases. Utilizing GPS for position and attitude determination, 
Kalman filters are used to determine all desired outputs. It is called by GOMONTE, 
which is basically a loop to perform the rendezvous for Monte Carlo simulations. 

RVSVLOS - This function determines which SVs in the constellation are in view. 

Inputs: 

x - current true state 

svpos - the row vector containing the positions of the GPS SVs, km (SV lx , 
SV ly , SV 1Z , SV 2x , SV 2y , SV 2z , ...) 

elevang - the desired elevation angle from the earth tangent, deg 
alt - the altitude of the rendezvous, km 
OUTPUTS: 
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svlos - the line of sight for all SVs, I = in sight, 0 = out of sight 


5A - This program calculates the selective availability residuals via the Braasch method. 
Inputs: 

blocklen - the length of time to create SA, sec 
numblock - the number of time histories to create 
userdt - the increment to output the SA in seconds, i.e., 60 = 1 min 
increments, 120 = 2 min increments 
Outputs: 

sares - matrix of SA residuals (timexblocks) 
sares2 - matrix of SA residuals (timexblocks) 
tsares - time index based on userdt 

STATNOIS - This function calculates the state noise for the Kalman filter, based on a 
specified percent unknown in the acceleration due to drag at the rendezvous altitude. 
Inputs: 

w - orbital angular velocity, rad/sec 
r - rendezvous altitude, km 
dt - measurement intervals, sec 
dragpct - the percent of acceleration due to drag to use 
Outputs: 

possn - position state noise, km 
velsn - velocity state noise , km/sec 

STMATRIX - This function calculates the chaser state transition matrix based on the 
linearized relative equations of motion. 

Inputs: 

t - measurement intervals, sec 
w - orbital angular velocity, rad/sec 
Outputs: 

stm - state transition matrix over the time interval t 

TRMOBSM - This function calculates the observation matrix for use in the Kalman 
filtering during the stationkeeping and terminal portions of the rendezvous. 
Inputs: 

h - the matrix of unit vectors to the SVs (SV,, SV 2 , SV v ... x x, y, z) 
Outputs: 

H - the observation matrix 

VOLTETRA - This function calculates the volume of the tetrahedron formed from four 
points. 

Inputs: 

a, b, c, d - the four points forming the tetrahedron 
Outputs: 

x - the associated volume 
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